using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

using RSCore.Helper.Collection;
using RSEngine.Math;

namespace RSEngine.Core
{

    sealed public class Mesh
    {

        #region Variables

        public const int F0DEL = 1 << 16;
        public const int F1DEL = 2 << 16;
        public const int SINGLEFACE = 3 << 16;
        public const int DELETED = 4 << 16;
        public const int SPLIT = 8 << 16;
        public const int BOUNDARY = 16 << 16;
        public const int LABELMASK = 65535;
        public const int BOUNDARY_RECTANGLE = 1;
        public const int BOUNDARY_CIRCLE = 2;
        public const int BOUNDARY_CHORD_RECTANGLE = 3;
        public const int BOUNDARY_CHORD_CIRCLE = 4;
        public const int BOUNDARY_PRESET = 5;
        public const int GEOMETRY_IMAGE = 1;
        public const int NORMAL_MAP = 2;

        public UnicityList<Vertex> Vertices = new UnicityList<Vertex>();
        public UnicityList<Edge> Edges = new UnicityList<Edge>();
        public UnicityList<Face> Faces = new UnicityList<Face>();

        #endregion

        #region Constructors

        public Mesh(Mesh M)
        {
            Copy(M);
        }

        public Mesh()
        {
        }

        #endregion

        public bool test(int i)
        {
            if (i == 0) return false;
            else return true;

            List<Vertex> l;
        }

        #region Copy

        public void Copy(Mesh M)
        {
            Destroy();
            for (int i = 0; i < M.Vertices.Count; i++)
            {
                M.Vertices[i].IndexInMesh = AddVertex(M.Vertices[i].Coordinate);
                Vertices[M.Vertices[i].IndexInMesh].TexCoordinate = (Vector3)M.Vertices[i].TexCoordinate.Clone();
                Vertices[M.Vertices[i].IndexInMesh].Normal = (Vector3)M.Vertices[i].Normal.Clone();
                Vertices[M.Vertices[i].IndexInMesh].n = (Vector3)M.Vertices[i].n.Clone();
            }
            for (int i = 0; i < M.Faces.Count; i++)
            {
                AddTriangle(M.Faces[i].Vertices[0].IndexInMesh,
                    M.Faces[i].Vertices[1].IndexInMesh,
                    M.Faces[i].Vertices[2].IndexInMesh);
            }
        }

        #endregion

        #region Clone

        public Mesh Clone()
        {
            Mesh m = new Mesh();
            m.Copy(this);
            return m;
        }

        #endregion

        #region Destroy

        public void Destroy()
        {
            Vertices.Clear();
            Edges.Clear();
            Faces.Clear();
        }

        #endregion

        //will always return vertices in sequence, unless vertices are deleted
        public int AddVertex(Vector3 v)
        {
            Vertex vert = new Vertex();
            vert.Coordinate = (Vector3)v.Clone();
            vert.n = new Vector3(0, 1, 0);
            vert.Normal = new Vector3(0, 1, 0);
            vert.TexCoordinate = new Vector3(0, 0, 0);
            return Vertices.AddWithoutCheck(vert);
        }

        public void SetTex(int v, Vector3 tex)
        {
            if (v < 0) return;
            if (v >= Vertices.Count) return;
            Vertices[v].TexCoordinate = (Vector3)tex.Clone();
        }

        public void SetVertex(int v, Vector3 vert)
        {
            if (v < 0) return;
            if (v >= Vertices.Count) return;
            Vertices[v].Coordinate = (Vector3)vert.Clone();
        }

        public void SetNormal(int v, Vector3 n)
        {
            if (v < 0) return;
            if (v >= Vertices.Count) return;
            Vertices[v].n = (Vector3)n.Clone();
            Vertices[v].Normal = (Vector3)n.Clone();
        }

        public Vertex GetVertex(int i)
        {
            if (i < 0) return null;
            if (i >= Vertices.Count) return null;
            return Vertices[i];
        }

        public void Recenter()
        {
            Vector3 min, max, c;
            if (Vertices.Count == 0)
                return;

            min = (Vector3)Vertices[0].Coordinate.Clone();
            max = (Vector3)Vertices[0].Coordinate.Clone();
            for (int i = 1; i < Vertices.Count; i++)
            {
                if (Vertices[i].Coordinate.x < min.x)
                    min.x = Vertices[i].Coordinate.x;
                if (Vertices[i].Coordinate.y < min.y)
                    min.y = Vertices[i].Coordinate.y;
                if (Vertices[i].Coordinate.z < min.z)
                    min.z = Vertices[i].Coordinate.z;
                if (Vertices[i].Coordinate.x > max.x)
                    max.x = Vertices[i].Coordinate.x;
                if (Vertices[i].Coordinate.y > max.y)
                    max.y = Vertices[i].Coordinate.y;
                if (Vertices[i].Coordinate.z > max.z)
                    max.z = Vertices[i].Coordinate.z;
            }

            c = (max + min) / 2.0;

            for (int i = 0; i < Vertices.Count; i++)
                Vertices[i].Coordinate -= c;
        }

        public void Rescale()
        {
            int i;
            Vector3 min, max;
            double maxw;
            if (Vertices.Count == 0) return;
            min = (Vector3)Vertices[0].Coordinate.Clone();
            max = (Vector3)Vertices[0].Coordinate.Clone();
            for (i = 1; i < Vertices.Count; i++)
            {
                if (Vertices[i].Coordinate.x < min.x)
                    min.x = Vertices[i].Coordinate.x;
                if (Vertices[i].Coordinate.y < min.y)
                    min.y = Vertices[i].Coordinate.y;
                if (Vertices[i].Coordinate.z < min.z)
                    min.z = Vertices[i].Coordinate.z;
                if (Vertices[i].Coordinate.x > max.x)
                    max.x = Vertices[i].Coordinate.x;
                if (Vertices[i].Coordinate.y > max.y)
                    max.y = Vertices[i].Coordinate.y;
                if (Vertices[i].Coordinate.z > max.z)
                    max.z = Vertices[i].Coordinate.z;
            }
            maxw = max.x - min.x;
            if (max.y - min.y > maxw) maxw = max.y - min.y;
            if (max.z - min.z > maxw) maxw = max.z - min.z;
            for (i = 0; i < Vertices.Count; i++)
            {
                Vertices[i].Coordinate /= maxw;
            }
        }

        public int AddTriangle(int v1, int v2, int v3)
        {
            int i;
            Vertex v;
            if ((v1 == v2) || (v1 == v3) || (v2 == v3)) return -1;
            Edge e1 = AddEdge(v1, v2);
            Edge e2 = AddEdge(v2, v3);
            Edge e3 = AddEdge(v3, v1);
            Face f = new Face();
            f.Vertices.Add(Vertices[v1]);
            f.Vertices.Add(Vertices[v2]);
            f.Vertices.Add(Vertices[v3]);
            f.Edges.Add(e1);
            f.Edges.Add(e2);
            f.Edges.Add(e3);
            //check for duplicate triangles
            //if the face exists, all vertices will have the face
            //   in their adjacency lists
            v = Vertices[v1];
            for (i = 0; i < v.Faces.Count; i++)
            {
                if (f.Equals(v.Faces[i]))
                {
                    return Faces.IndexOf(v.Faces[i]);
                }
            }
            Vertices[v1].Faces.Add(f);
            Vertices[v2].Faces.Add(f);
            Vertices[v3].Faces.Add(f);
            e1.Faces.Add(f);
            e2.Faces.Add(f);
            e3.Faces.Add(f);
            return Faces.AddWithoutCheck(f);
        }

        //take a triangle already in the mesh and read it
        public void AddTriangle(int v1, int v2, int v3, Face f)
        {
            int i;
            Vertex v;
            Edge e1 = AddEdge(v1, v2);
            Edge e2 = AddEdge(v2, v3);
            Edge e3 = AddEdge(v3, v1);
            f.Vertices.Add(Vertices[v1]);
            f.Vertices.Add(Vertices[v2]);
            f.Vertices.Add(Vertices[v3]);
            f.Edges.Add(e1);
            f.Edges.Add(e2);
            f.Edges.Add(e3);
            //check for duplicate triangles
            //if the face exists, all vertices will have the face
            //   in their adjacency lists
            v = Vertices[v1];
            for (i = 0; i < v.Faces.Count; i++)
            {
                if (f.Equals(v.Faces[i]))
                {
                    return;
                }
            }
            Vertices[v1].Faces.Add(f);
            Vertices[v2].Faces.Add(f);
            Vertices[v3].Faces.Add(f);
            e1.Faces.Add(f);
            e2.Faces.Add(f);
            e3.Faces.Add(f);
        }

        public void ReverseNormals()
        {
            int i;

            for (i = 0; i < Vertices.Count; i++)
            {
                Vertices[i].Normal = -Vertices[i].Normal;
            }
            for (i = 0; i < Faces.Count; i++)
            {
                Faces[i].Normal = -Faces[i].Normal;
            }
        }

        public int Genus()
        {
            int i, j;
            int N, F, E, G, B;
            bool loop = true;
            //Euler-Poincare formula: N+F-E=2-2G-B
            //N=no. vertices
            //F=no. faces
            //E=no. edges
            //G=genus
            //B=no. of boundary loops
            //thus G=(2-B-N-F+E)/2;

            N = Vertices.Count;
            F = Faces.Count;
            E = Edges.Count;
            for (i = 0; i < E; i++)
            {
                Edges[i].label = 0;
            }
            B = 0;
            //visit boundaries
            while (loop)
            {
                loop = false;
                for (i = 0; i < E; i++)
                {
                    if ((Edges[i].label == 0) && (Edges[i].IsBoundary))
                    {
                        loop = true;
                        B++;
                        Edge e = Edges[i];
                        while (e != null)
                        {
                            e.label = 1;
                            Edge e2 = null;
                            for (j = 0; j < e.Vertices[0].Edges.Count; j++)
                            {
                                if ((e.Vertices[0].Edges[j].label == 0) && (e.Vertices[0].Edges[j].IsBoundary))
                                {
                                    e2 = e.Vertices[0].Edges[j];
                                }
                            }
                            if (e2 == null)
                                for (j = 0; j < e.Vertices[1].Edges.Count; j++)
                                {
                                    if ((e.Vertices[1].Edges[j].label == 0) && (e.Vertices[1].Edges[j].IsBoundary))
                                    {
                                        e2 = e.Vertices[1].Edges[j];
                                    }
                                }
                            e = e2;
                        }
                    }
                    else Edges[i].label = 1;
                }
            }

            G = (2 - B - N - F + E);
            return G / 2;
        }

        public void ComputeFaceNormals()
        {
            int i;
            for (i = 0; i < Faces.Count; i++)
                Faces[i].UpdateNormal();
        }

        public void ComputeVertexNormals()
        {
            int i;

            ComputeFaceNormals();
            for (i = 0; i < Vertices.Count; i++)
                Vertices[i].UpdateNormal();
        }

        public bool Connected()
        {
            if (!ConnectedFaces())
                return false;
            if (!ConnectedEdges())
                return false;

            return true;
        }

        public bool ConnectedFaces()
        {
            int i;

            if (Faces.Count == 0) return true;
            for (i = 0; i < Faces.Count; i++)
                Faces[i].label = 0;

            CheckConnected(Faces[0]);
            for (i = 0; i < Faces.Count; i++)
            {
                if (Faces[i].label == 0) return false;
            }
            return true;
        }

        public bool ConnectedEdges()
        {
            int i;

            if (Edges.Count == 0) return true;

            for (i = 0; i < Edges.Count; i++)
            {
                Edges[i].label = 0;
            }

            for (i = 0; i < Vertices.Count; i++)
            {
                Vertices[i].label = 0;
            }

            CheckConnected(Edges[0]);
            for (i = 0; i < Edges.Count; i++)
            {
                if (Edges[i].label == 0) return false;
            }
            for (i = 0; i < Vertices.Count; i++)
            {
                if (Vertices[i].label == 0) return false;
            }
            return true;
        }

        //only 2-manifold meshes can be cut and parameterised.
        public bool Manifold()
        {
            int i;
            for (i = 0; i < Edges.Count; i++)
            {
                if (Edges[i].Faces.Count > 2) return false;
            }
            return true;
        }

        // Check that no vertex is connected to two or more
        // triangle fans, if so, split the fans
        // This happens if the mesh touches itself at exactly
        // one point, a vertex (i.e. not an edge or face). 
        public void CheckFans()
        {
            int i, j, k;
            bool recheck = false;
            for (i = 0; i < Vertices.Count; i++)
            {
                Vertex v = Vertices[i];
                for (j = 0; j < v.Faces.Count; j++)
                {
                    v.Faces[j].label = 0;
                }
                if (v.Faces.Count == 0) continue;
                v.Faces[0].label = 1;
                bool changed = true;
                while (changed)
                {
                    changed = false;
                    for (j = 0; j < v.Faces.Count; j++)
                    {
                        for (k = j + 1; k < v.Faces.Count; k++)
                        {
                            if ((v.Faces[j].label == 0) || (v.Faces[k].label == 0))
                            {
                                Edge e = v.Faces[j].FindAdjacent(v.Faces[k]);
                                if (e != null)
                                {
                                    if (v.Faces[j].label > 0)
                                    {
                                        v.Faces[k].label = v.Faces[j].label;
                                        changed = true;
                                    }
                                    else if (v.Faces[k].label > 0)
                                    {
                                        v.Faces[j].label = v.Faces[k].label;
                                        changed = true;
                                    }
                                }
                            }
                        }
                    }
                }

                bool touch = false;
                for (j = 0; j < v.Faces.Count; j++)
                {
                    if (v.Faces[j].label == 0)
                    {
                        touch = true;
                    }
                }
                if (touch == true)
                {
                    Console.Write("Found touching vertex\n");
                    Vertex vnew = new Vertex();
                    vnew.Distance = v.Distance;
                    vnew.label = v.label;
                    vnew.n = (Vector3)v.n.Clone();
                    vnew.Normal = (Vector3)v.Normal.Clone();
                    vnew.TexCoordinate = (Vector3)v.TexCoordinate.Clone();
                    vnew.Coordinate = (Vector3)v.Coordinate.Clone();
                    UnicityList<Face> vfaces = v.Faces.Clone();
                    Vertices.AddWithoutCheck(vnew);
                    recheck = true;
                    for (j = 0; j < vfaces.Count; j++)
                    {
                        if (vfaces[j].label == 1)
                        {
                            vfaces[j].Replace(v, vnew);
                            for (k = 0; k < vfaces[j].Edges.Count; k++)
                            {
                                vfaces[j].Edges[k].Replace(v, vnew);
                            }
                        }
                    }
                }
            }

            if (recheck) CheckFans();
        }

        // Cut mesh to be equivalent to a disc
        // only works if the mesh does not
        // "cross" itself. That is, inside and outside is always clear.
        public void DiscCut()
        {
            int i, j;
            UnicityList<Edge> cutpath = new UnicityList<Edge>();
            UnicityList<Edge> candidates = new UnicityList<Edge>();
            Edge selected;
            double min;
            int valence, x;
            Vertex v;
            bool changed;
            bool boundary = false;

            //mark boundary edges
            for (i = 0; i < Edges.Count; i++)
            {
                Edges[i].label = 0;
                //mark boundaries, marked as split as well, so that we don't try to split them
                if (Edges[i].Faces.Count < 2)
                {
                    Edges[i].label |= F1DEL | BOUNDARY | SPLIT;
                    boundary = true;
                }
            }

            //no faces are deleted
            for (i = 0; i < Faces.Count; i++)
            {
                Faces[i].label = 0;
            }

            //for computation of geodesic distance
            for (i = 0; i < Vertices.Count; i++)
            {
                Vertices[i].Distance = 1e15;
            }
            //remove a seed triangle
            int fs = 0;
            DiscDeleteFace(Faces[fs], candidates);
            //update geodesic distance for each 
            // of the seed triangle vertices
            UnicityList<Vertex> border = new UnicityList<Vertex>();
            border.Add(Faces[fs].Vertices[0]);
            border.Add(Faces[fs].Vertices[1]);
            border.Add(Faces[fs].Vertices[2]);
            Console.Write("Distance...\n");
            ComputeGeodesicDistance(border);

            //successively remove edges that are only adjacent 
            // to one triangle, and remove the edge and triangle
            Console.Write("Removing triangles...\n");
            selected = candidates[0];
            while (selected != null)
            {
                selected = null;
                min = 1e16;
                min = -1.0;
                for (i = 0; i < candidates.Count; i++)
                {
                    if ((!(test(candidates[i].label & BOUNDARY))) &&
                        ((candidates[i].label & SINGLEFACE) > 0))
                    {
                        if (!((test(candidates[i].label & F0DEL)) && (test(candidates[i].label & F1DEL))))
                        {
                            if (candidates[i].Vertices[0].Distance > min)
                            {
                                selected = candidates[i];
                                min = selected.Vertices[0].Distance;
                            }
                            if (candidates[i].Vertices[1].Distance > min)
                            {
                                selected = candidates[i];
                                min = selected.Vertices[1].Distance;
                            }
                        }
                    }
                }
                if (selected != null)
                {
                    //one of the two must be deleted,
                    // one is deleted already
                    //  DiscDeleteFace checks if the face is already deleted
                    // The corresponding edge is also deleted.
                    candidates.Remove(selected);
                    DiscDeleteEdge(selected);
                    DiscDeleteFace(selected.Faces[0], candidates);
                    DiscDeleteFace(selected.Faces[1], candidates);
                }
            }

            //the cut path can now be built
            for (i = 0; i < Edges.Count; i++)
            {
                if ((!(test(Edges[i].label & DELETED))) &&
                    (!(test(Edges[i].label & BOUNDARY))))
                    cutpath.Add(Edges[i]);
            }

            Console.Write("Pruning...{0}\n", cutpath.Count);
            //prune vertices that
            // are connected to only one edge, delete the edge also.
            // Repeat until only loops remain
            changed = true;
            while (changed)
            {
                changed = false;
                for (i = cutpath.Count - 1; i >= 0; i--)
                {
                    valence = 0;
                    v = cutpath[i].Vertices[0];
                    for (j = 0; j < v.Edges.Count; j++)
                    {
                        if (!(test(v.Edges[j].label & DELETED))) valence++;
                    }
                    //store the valence in the label
                    v.label = valence;
                    //if the valence is one, prune the edge
                    if (valence == 1)
                    {
                        cutpath[i].label |= DELETED;
                        cutpath.Remove(i);
                        changed = true;
                        continue;
                    }
                    //Recompute valence by disregarding boundary edges, since they are not part of the cut
                    //  homever, if a boundary is found, then the valence should be increased by one, since
                    //  a boundary can be regarded as a cut that has already been made.
                    //  A vertex can be joined to a maximum of two boundary edges, any more and the graph is
                    //	not connected nicely. Should this be detected and cut elsewhere? 
                    valence = 0;
                    x = 0;
                    for (j = 0; j < v.Edges.Count; j++)
                    {
                        if (test(v.Edges[j].label & BOUNDARY)) x = 1;
                        if ((!(test(v.Edges[j].label & DELETED))) && (!(test(v.Edges[j].label & BOUNDARY)))) valence++;
                    }
                    //store the valence in the label
                    v.label = valence + x;

                    valence = 0;
                    v = cutpath[i].Vertices[1];
                    for (j = 0; j < v.Edges.Count; j++)
                    {
                        if (!(test(v.Edges[j].label & DELETED))) valence++;
                    }
                    //store the valence in the label
                    v.label = valence;
                    //if the valence is one, prune the edge
                    if (valence == 1)
                    {
                        cutpath[i].label |= DELETED;
                        cutpath.Remove(i);
                        changed = true;
                        continue;
                    }
                    //recompute valence by disregarding boundary edges, since they are not part of the cut
                    valence = 0;
                    x = 0;
                    for (j = 0; j < v.Edges.Count; j++)
                    {
                        if (test(v.Edges[j].label & BOUNDARY)) x = 1;
                        if ((!(test(v.Edges[j].label & DELETED))) && (!(test(v.Edges[j].label & BOUNDARY)))) valence++;
                    }
                    //store the valence in the label
                    v.label = valence + x;
                }
            }

            //now try to remove serrations in the path
            // we do this by checking the triangles adjacent to
            // each edge in the path, if another edge on the
            // path is part of the triangle, then we check 
            // the third edge. If the length of the third edge
            // is shorter than the first two, then we replace 
            // the two edges by the third edge.
            Console.Write("Smoothing cut...\n");
            while (SmoothCut(cutpath)) ;
            Random rnd = new Random();

            //Test for genus 0 surface
            if ((cutpath.Count == 0) && (!boundary))
            {
                //If we get a genus 0 surface, we should rather use
                // an alternative parameterisation and avoid the
                //	disc cut

                //create a random connected path to cut along	
                UnicityList<Vertex> cutvert = new UnicityList<Vertex>();
                i = rnd.Next() % Edges.Count;
                cutpath.Add(Edges[i]);
                min = Edges.Count;
                //min=sqrt(min)*2;
                min = 8;
                v = Edges[i].Vertices[1];
                cutvert.Add(Edges[i].Vertices[0]);
                cutvert.Add(Edges[i].Vertices[1]);
                //remark edges
                for (i = 0; i < Edges.Count; i++)
                {
                    Edges[i].label = DELETED;
                }
                for (i = 0; i < min; i++)
                {
                    selected = v.Edges[rnd.Next() % (v.Edges.Count)];
                    while (cutpath.Contains(selected))
                    {
                        selected = v.Edges[rnd.Next() % (v.Edges.Count)];
                    }
                    Vertex nv;
                    if (selected.Vertices[0] == v) nv = selected.Vertices[1];
                    else nv = selected.Vertices[0];
                    if (!cutvert.Contains(nv))
                    {
                        cutvert.Add(nv);
                        cutpath.Add(selected);
                        v = nv;
                    }
                }
                //remark edges
                for (i = 0; i < cutpath.Count; i++)
                {
                    cutpath[i].label = 0;
                }
                //compute valences for path
                for (i = 0; i < cutvert.Count; i++)
                {
                    valence = 0;
                    v = cutvert[i];
                    for (j = 0; j < v.Edges.Count; j++)
                    {
                        if (!(test(v.Edges[j].label & DELETED))) valence++;
                    }
                    //store the valence in the label
                    v.label = valence;
                }
                //while (SmoothCut(cutpath));
            }
            if ((cutpath.Count == 0) && (!boundary))
            {
                //failed to cut the surface!
                //we should hopefully never get here
            }
            Console.Write("Performing cut...{0}\n", cutpath.Count);
            CutMesh(cutpath);
            Console.Write("Cut done.\n");
        }

        // Cut for geometry image using the technique described in Gu's thesis 
        // on geometry images
        // C# doesn't support default parameters, so two 
        // methods are required
        public void GuCut()
        {
            GuCut(false);
        }

        public void GuCut(bool flipedge)
        {
            UnicityList<Edge> eboundary = new UnicityList<Edge>();
            UnicityList<Vertex> vboundary = new UnicityList<Vertex>();
            Edge e, ce;
            Vertex v, vp;
            Face f;
            int i, j, c;
            double maxstretch, stretch, Eold, Enew;
            UnicityList<Edge> cutpath = new UnicityList<Edge>();
            Mesh opt = null;
            int iter = 1000;

            //first cut into a disc
            Console.Write("Cutting\n");
            DiscCut();

            //then compute Floater parameterisation
            Console.Write("Floater param to circle\n");
            ParamFloater(BOUNDARY_CIRCLE);
            Enew = AllStretch();
            Eold = Enew + 1.0;
            Console.Write("Stretch before: {0}\n", Enew);
            //refine until stretch increases
            while ((Eold > Enew) && (--iter > 0))
            {
                //save curent mesh
                opt = this.Clone();
                if (!flipedge)
                {
                    //If any border triangle is degenerate in parameter space,
                    // that is, three vertices are on the border, split all edges not on the border.
                    // This creates  4 new triangles (usually), and 2 are discarded 
                    cutpath.Clear();
                    for (i = 0; i < Faces.Count; i++)
                    {
                        c = 0;
                        for (j = 0; j < 3; j++)
                            if (Faces[i].Vertices[j].IsBoundary)
                                c++;
                        if (c == 3)
                        {
                            for (j = 0; j < 3; j++)
                                if (!(Faces[i].Edges[j].IsBoundary))
                                {
                                    cutpath.Add(Faces[i].Edges[j]);
                                }
                        }
                    }
                    if (cutpath.Count > 0)
                    {
                        for (i = 0; i < cutpath.Count; i++)
                        {
                            SplitEdge(cutpath[i]);
                        }
                        //reparameterise
                        Console.Write("Parameterise...\n");
                        ParamFloater(BOUNDARY_CIRCLE);
                        //check if stretch is minimised
                        Eold = Enew;
                        Enew = AllStretch();
                        Eold = Enew + 1.0;
                        Console.Write("Stretch new: {0}\n", Enew);
                        continue;
                    }
                }
                else
                {
                    //flip the offending edge instead of introducing new vertices
                    for (i = 0; i < Vertices.Count; i++)
                    {
                        Vertices[i].IndexInMesh = i;
                    }
                    for (i = 0; i < Faces.Count; i++)
                    {
                        c = 0;
                        for (j = 0; j < 3; j++)
                            if (Faces[i].Vertices[j].IsBoundary)
                                c++;
                        if (c == 3)
                        {
                            Edge et = null;
                            for (j = 0; j < 3; j++)
                            {
                                if (!Faces[i].Edges[j].IsBoundary)
                                {
                                    et = Faces[i].Edges[j];
                                }
                            }
                            if (et != null)
                            {
                                Vertex v1, v2, v3 = null, v4 = null;
                                v1 = et.Vertices[0];
                                v2 = et.Vertices[1];
                                for (j = 0; j < 3; j++)
                                {
                                    if ((et.Faces[0].Vertices[j] != v1) &&
                                        (et.Faces[0].Vertices[j] != v2))
                                        v3 = et.Faces[0].Vertices[j];
                                    if ((et.Faces[1].Vertices[j] != v1) &&
                                        (et.Faces[1].Vertices[j] != v2))
                                        v4 = et.Faces[1].Vertices[j];
                                }
                                RemoveFace(et.Faces[0]);
                                RemoveFace(et.Faces[1]);
                                Edges.Remove(et);
                                AddTriangle(v1.IndexInMesh, v3.IndexInMesh, v4.IndexInMesh);
                                AddTriangle(v2.IndexInMesh, v3.IndexInMesh, v4.IndexInMesh);
                            }
                        }
                    }

                    //reparameterise
                    Console.Write("Parameterise...\n");
                    ParamFloater(BOUNDARY_CIRCLE);
                    //check if stretch is minimised
                    Eold = Enew;
                    Enew = AllStretch();
                    Eold = Enew + 1.0;
                    Console.Write("Stretch new: {0}\n", Enew);
                    continue;
                }
                maxstretch = -1;
                f = null;
                for (i = 0; i < Faces.Count; i++)
                {
                    //don't consider triangles on the boundary
                    if (Faces[i].Edges[0].IsBoundary) continue;
                    if (Faces[i].Edges[1].IsBoundary) continue;
                    if (Faces[i].Edges[2].IsBoundary) continue;
                    stretch = Faces[i].Stretch();
                    if (maxstretch < stretch)
                    {
                        maxstretch = stretch;
                        f = Faces[i];
                    }
                }
                if (f == null) Console.Write("No maximal stretch!\n");
                //Check for high stretch.
                //if found, join vertex of high stretch to current border 
                //  using shortest path	
                eboundary.Clear();
                vboundary.Clear();
                ComputeBoundary(eboundary, vboundary);

                Console.Write("Calculating distances...\n");
                ComputeGeodesicDistance(vboundary);

                v = f.Vertices[0];
                cutpath.Clear();
                //follow path from vertex to boundary
                Console.Write("Building cutpath...\n");
                while (v.Next != null)
                {
                    vp = v.Next;
                    ce = new Edge();
                    ce.Vertices[0] = v;
                    ce.Vertices[1] = vp;
                    //find matching edge
                    e = v.Edges.Match(ce);
                    if (e == null) Console.Write("Bad edge! v1={0} v2={1}\n", v, vp);
                    //add to cutpath
                    cutpath.Add(e);
                    if (e.IsBoundary) Console.Write("Boundary edge!\n");
                    //advance to next node
                    v = vp;
                }

                Console.Write("Cutting...\n");
                //cut the mesh
                CutMesh(cutpath);

                //reparameterise
                Console.Write("Parameterise...\n");
                ParamFloater(BOUNDARY_CIRCLE);
                //check if stretch is minimised
                Eold = Enew;
                Enew = AllStretch();
                Console.Write("Stretch new: {0}\n", Enew);
            }
            //revert to previous mesh (stretch was better)
            Copy(opt);
            Enew = AllStretch();
            Console.Write("Stretch after: {0}\n", Enew);
        }

        public void SplitEdge(Edge e)
        {
            Face f;
            Vector3 pos;
            int v1, v2, v3, vpos, i, j;

            pos = (e.Vertices[0].Coordinate + e.Vertices[1].Coordinate) / 2.0;
            vpos = AddVertex(pos);

            //remove edge
            e.Vertices[0].Edges.Remove(e);
            e.Vertices[1].Edges.Remove(e);
            Edges.Remove(e);

            for (i = 0; i < 2; i++)
            {
                f = e.Faces[i];
                if (f != null)
                {
                    for (j = 0; j < 3; j++)
                    {
                        f.Vertices[j].Faces.Remove(f);
                        f.Edges[j].Faces.Remove(f);
                    }
                    for (j = 0; j < 3; j++)
                    {
                        if ((f.Vertices[j] != e.Vertices[0]) && (f.Vertices[j] != e.Vertices[1]))
                        {
                            Vertex tmp = f.Vertices[0];
                            f.Vertices[0] = f.Vertices[j];
                            f.Vertices[j] = tmp;
                        }
                    }
                    v1 = Vertices.IndexOf(f.Vertices[0]);
                    v2 = Vertices.IndexOf(f.Vertices[1]);
                    v3 = Vertices.IndexOf(f.Vertices[2]);

                    Faces.Remove(f);
                    AddTriangle(v1, v2, vpos);
                    AddTriangle(v1, vpos, v3);
                }
            }

        }

        public void CutMesh(UnicityList<Edge> cutpath)
        {
            int i, j, k, valence, x;
            UnicityList<Vertex> cutvert = new UnicityList<Vertex>();
            Vertex[] newvert;
            int count;
            Vertex v;
            Edge adj;
            bool changed;
            Edge e;

            Console.Write("Cutting {0}/{1} edges\n", cutpath.Count, Edges.Count);
            //label edges as if they are not in the cutpath
            for (i = 0; i < Edges.Count; i++)
            {
                if (Edges[i].IsBoundary) Edges[i].label = BOUNDARY | SPLIT;
                else Edges[i].label = DELETED;
            }

            //label edges in cutpath as not deleted
            for (i = 0; i < cutpath.Count; i++)
            {
                //label edge as not deleted
                if (!cutpath[i].IsBoundary)
                    cutpath[i].label = 0;
                else
                {
                    Console.Write("Boundary in cutpath!\n");
                    cutpath.Remove(i);
                    i--;
                }
            }

            //label vertices according to valence (undeleted edges)
            for (i = 0; i < cutpath.Count; i++)
            {
                valence = 0;
                v = cutpath[i].Vertices[0];
                x = 0;
                for (j = 0; j < v.Edges.Count; j++)
                {
                    if (test(v.Edges[j].label & BOUNDARY)) x = 1;
                    if ((!(test(v.Edges[j].label & DELETED))) && (!(test(v.Edges[j].label & BOUNDARY)))) valence++;
                }
                //store the valence in the label
                v.label = valence + x;

                valence = 0;
                v = cutpath[i].Vertices[1];
                x = 0;
                for (j = 0; j < v.Edges.Count; j++)
                {
                    if (test(v.Edges[j].label & BOUNDARY)) x = 1;
                    if ((!(test(v.Edges[j].label & DELETED))) && (!(test(v.Edges[j].label & BOUNDARY)))) valence++;
                }
                //store the valence in the label
                v.label = valence + x;
            }

            //find the vertices in the cut
            for (i = 0; i < cutpath.Count; i++)
            {
                //check valence
                //for genus 0 surfaces, the cut will not necessarily be a loop
                if (cutpath[i].Vertices[0].label >= 2)
                    cutvert.Add(cutpath[i].Vertices[0]);
                if (cutpath[i].Vertices[1].label >= 2)
                    cutvert.Add(cutpath[i].Vertices[1]);
            }

            for (i = 0; i < cutvert.Count; i++)
            {
                count = 0;
                v = cutvert[i];
                //create a new vertex for each edge
                // entering this vertex (each edge on a cutpath)
                newvert = new Vertex[v.label];
                newvert[0] = v;
                for (j = 1; j < v.label; j++)
                {
                    newvert[j] = new Vertex();
                    newvert[j].Coordinate = (Vector3)v.Coordinate.Clone();
                    newvert[j].TexCoordinate = (Vector3)v.TexCoordinate.Clone();
                    newvert[j].Normal = (Vector3)v.Normal.Clone();
                    newvert[j].n = (Vector3)v.n.Clone();
                    newvert[j].Distance = v.Distance;
                    Vertices.Add(newvert[j]);
                }
                //reset labels for cut
                for (j = 0; j < v.Faces.Count; j++)
                {
                    v.Faces[j].label = -1;
                }

                changed = true;
                //label the faces
                while (changed)
                {
                    changed = false;
                    //first label a triangle next to a cut edge
                    for (j = 0; (!changed) && (j < v.Faces.Count); j++)
                    {
                        if (v.Faces[j].label < 0)
                            for (k = 0; (!changed) && (k < 3); k++)
                            {
                                if ((!(test(v.Faces[j].Edges[k].label & DELETED)))
                                    && (!(test(v.Faces[j].Edges[k].label & BOUNDARY)))
                                    //	&&(!(v.Faces[j].edges[k].label&SPLIT)) 
                                )
                                {
                                    v.Faces[j].label = count;
                                    changed = true;
                                    count++;
                                }
                            }
                    }

                    if (changed == false) break;
                    //now propogate the information to the surrounding faces
                    while (changed)
                    {
                        changed = false;
                        for (j = 0; j < v.Faces.Count; j++)
                        {
                            for (k = j + 1; k < v.Faces.Count; k++)
                            {
                                if ((v.Faces[j].label < 0) || (v.Faces[k].label < 0))
                                {
                                    adj = v.Faces[j].FindAdjacent(v.Faces[k]);
                                    if (adj != null)
                                    {
                                        if (test(adj.label & DELETED))
                                        {
                                            //faces should have same label
                                            if (v.Faces[j].label >= 0)
                                            {
                                                v.Faces[k].label = v.Faces[j].label;
                                                changed = true;
                                            }
                                            else if (v.Faces[k].label >= 0)
                                            {
                                                v.Faces[j].label = v.Faces[k].label;
                                                changed = true;
                                            }
                                        }
                                        else
                                        {
                                            //faces should have different labels
                                            //not sure what yet though
                                        }
                                    }
                                }
                            }
                        }
                    }
                    //go to next iteration, to choose a new cut face to label
                    changed = true;
                }

                for (j = 0; j < v.Faces.Count; j++)
                {
                    if (v.Faces[j].label < 0)
                    {
                        Console.Write("Critical error, unlabeled face\n");
                    }
                    if (v.Faces[j].label >= v.label)
                    {
                        Console.Write("Critical error, incorrectly labeled face v.label={0} face={1}\n", v.label, v.Faces[j].label);
                    }
                }

                for (j = 0; j < v.Edges.Count; j++)
                {
                    if (!(test(v.Edges[j].label & DELETED)))
                    {
                        if (v.Edges[j].Faces.Count >= 2)
                        {
                            if (v.Edges[j].Faces[0].label == v.Edges[j].Faces[1].label)
                            {
                                Console.Write("Critical error, faces have same label but shouldn't\n");
                                Console.Write("label={0}\n", v.Edges[j].Faces[0].label);
                            }
                        }
                    }
                    else
                    {
                        if (v.Edges[j].Faces.Count >= 2)
                        {
                            if (v.Edges[j].Faces[0].label != v.Edges[j].Faces[1].label)
                            {
                                Console.Write("Critical error, faces should have the same label\n");
                            }
                        }
                    }
                }

                //now, update the edges to use the correct vertex
                //if the edge was a cut edge, split the edge and use the
                //   face labels to determine which edge belongs
                //   to which face. These new edges are no longer cut edges.

                //Console.Write("Splitting edges\n");
                // go backwards, because we will be removing edges while working	
                for (j = v.Edges.Count - 1; j >= 0; j--)
                {
                    if ((test(v.Edges[j].label & DELETED)) || (test(v.Edges[j].label & SPLIT)))
                    {
                        //edge not part of cutpath: just change the vertex
                        //Console.Write("Splitting non-cutpath edge\n");	
                        e = v.Edges[j];
                        if (e.Faces[0] != null)
                        {
                            if (e.Vertices[0] == v)
                            {
                                v.Edges.Remove(e);
                                e.Vertices[0] = newvert[e.Faces[0].label];
                                e.Vertices[0].Edges.Add(e);
                            }
                            else
                                if (e.Vertices[1] == v)
                                {
                                    v.Edges.Remove(e);
                                    e.Vertices[1] = newvert[e.Faces[0].label];
                                    e.Vertices[1].Edges.Add(e);
                                }
                        }
                        else
                        {
                            if (e.Vertices[0] == v)
                            {
                                v.Edges.Remove(e);
                                e.Vertices[0] = newvert[e.Faces[1].label];
                                e.Vertices[0].Edges.Add(e);
                            }
                            else
                                if (e.Vertices[1] == v)
                                {
                                    v.Edges.Remove(e);
                                    e.Vertices[1] = newvert[e.Faces[1].label];
                                    e.Vertices[1].Edges.Add(e);
                                }
                        }
                        //Console.Write("Split done\n");	
                    }
                    else
                    {
                        //change the existing edge to use the correct vertex
                        //Console.Write("Splitting cutpath edge face0=%x face1=%x\n", v.edges[j].faces[0], v.edges[j].faces[1]);	
                        Edge newedge;
                        newedge = v.Edges[j].Clone();
                        e = v.Edges[j];
                        if (e.Vertices[0] == v)
                        {
                            v.Edges.Remove(e);
                            e.Vertices[0] = newvert[e.Faces[0].label];
                            e.Vertices[0].Edges.Add(e);
                        }
                        else
                            if (e.Vertices[1] == v)
                            {
                                v.Edges.Remove(e);
                                e.Vertices[1] = newvert[e.Faces[0].label];
                                e.Vertices[1].Edges.Add(e);
                            }
                        //Console.Write("Actual split\n");	
                        //edge part of cutpath: split edge
                        //Console.Write("Replacing ends\n");
                        if (newedge.Vertices[0] == v)
                        {
                            newedge.Vertices[0] = newvert[newedge.Faces[1].label];
                        }
                        else
                            if (newedge.Vertices[1] == v)
                            {
                                newedge.Vertices[1] = newvert[newedge.Faces[1].label];
                            }
                        //Console.Write("Adding edge to vertices\n");	
                        newedge.Vertices[0].Edges.Add(newedge);
                        newedge.Vertices[1].Edges.Add(newedge);
                        //Console.Write("Selecting faces\n");	
                        //newedge.faces[0]=null;
                        newedge.Faces.Remove(0);
                        //don't have to set the face to null, since replace
                        // will do it for us
                        //e.faces[1]=null;
                        //Console.Write("Replacing edges in triangle %x\n", newedge.faces[1]);	
                        newedge.Faces[0].Replace(e, newedge);
                        newedge.label |= SPLIT;
                        e.label |= SPLIT;
                        //Console.Write("Adding edge to set of edges\n");	
                        Edges.Add(newedge);
                        //Console.Write("Split done\n");	
                    }
                }

                //now update the triangles
                for (j = v.Faces.Count - 1; j >= 0; j--)
                {
                    v.Faces[j].Replace(v, newvert[v.Faces[j].label]);
                }
            }
        }

        public bool SmoothCut(UnicityList<Edge> cutpath)
        {
            int i, j, k;
            Edge other;
            Face tri;
            Face f = new Face();
            int[] tag = new int[3] { 0, 0, 0 };
            Edge e1;
            Edge e2;
            Vertex v;
            int valence, x;

            for (i = 0; i < cutpath.Count; i++)
            {
                //check for later connecting edges
                for (j = i + 1; j < cutpath.Count; j++)
                {
                    //test if they share an endpoint, but there must be only two edges 
                    // sharing this vertex, thus the valence (stored in the label) must be 2
                    tri = null;
                    if (cutpath[i].Vertices[0].label == 2)
                    {
                        if (cutpath[j].Vertices[0] == cutpath[i].Vertices[0])
                        {
                            f.Vertices[0] = cutpath[i].Vertices[0];
                            f.Vertices[1] = cutpath[i].Vertices[1];
                            f.Vertices[2] = cutpath[j].Vertices[1];
                            tri = Faces.Match(f);
                        }
                        if (cutpath[j].Vertices[1] == cutpath[i].Vertices[0])
                        {
                            f.Vertices[0] = cutpath[i].Vertices[0];
                            f.Vertices[1] = cutpath[i].Vertices[1];
                            f.Vertices[2] = cutpath[j].Vertices[0];
                            tri = Faces.Match(f);
                        }
                    }
                    if (cutpath[i].Vertices[1].label == 2)
                    {
                        if (cutpath[j].Vertices[0] == cutpath[i].Vertices[1])
                        {
                            f.Vertices[0] = cutpath[i].Vertices[0];
                            f.Vertices[1] = cutpath[i].Vertices[1];
                            f.Vertices[2] = cutpath[j].Vertices[1];
                            tri = Faces.Match(f);
                        }
                        if (cutpath[j].Vertices[1] == cutpath[i].Vertices[1])
                        {
                            f.Vertices[0] = cutpath[i].Vertices[0];
                            f.Vertices[1] = cutpath[i].Vertices[1];
                            f.Vertices[2] = cutpath[j].Vertices[0];
                            tri = Faces.Match(f);
                        }
                    }
                    //if we find such a triangle, find the other edge
                    //  and see if it is shorter
                    if (tri != null)
                    {
                        for (k = 0; k < 3; k++) tag[k] = 0;
                        for (k = 0; k < 3; k++)
                        {
                            if (tri.Edges[k] == cutpath[j]) tag[k] = 1;
                            if (tri.Edges[k] == cutpath[i]) tag[k] = 1;
                        }
                        other = null;
                        for (k = 0; k < 3; k++)
                        {
                            if (tag[k] == 0) other = tri.Edges[k];
                        }
                        if (other.Length < cutpath[i].Length + cutpath[j].Length)
                        {
                            //the edge is shorter, replace the other two by this edge
                            e1 = cutpath[i];
                            e2 = cutpath[j];
                            cutpath.Add(other);
                            cutpath.Remove(i);
                            cutpath.Remove(j);
                            e1.label |= DELETED;
                            e2.label |= DELETED;
                            other.label &= ~DELETED;

                            //Recompute valence by disregarding boundary edges, since they are not part of the cut
                            //  homever, if a boundary is found, then the valence should be increased by one, since
                            //  a boundary can be regarded as a cut that has already been made.
                            //  A vertex can be joined to a maximum of two boundary edges, any more and the graph is
                            //	not connected nicely. Should this be detected and cut elsewhere? 
                            v = other.Vertices[0];
                            valence = 0;
                            x = 0;
                            for (j = 0; j < v.Edges.Count; j++)
                            {
                                if (test(v.Edges[j].label & BOUNDARY)) x = 1;
                                if ((!(test(v.Edges[j].label & DELETED))) && (!(test(v.Edges[j].label & BOUNDARY)))) valence++;
                            }
                            //store the valence in the label
                            v.label = valence + x;

                            v = other.Vertices[1];
                            valence = 0;
                            x = 0;
                            for (j = 0; j < v.Edges.Count; j++)
                            {
                                if (test(v.Edges[j].label & BOUNDARY)) x = 1;
                                if ((!(test(v.Edges[j].label & DELETED))) && (!(test(v.Edges[j].label & BOUNDARY)))) valence++;
                            }
                            //store the valence in the label
                            v.label = valence + x;

                            return true;
                        }
                    }
                }
            }
            return false;
        }

        // calculate the geodesic distance, the shortest path
        // on the surface, between this vertex and all other vertices
        // uses Dijkstra's algorithm.
        // The labels of vertices are set to point to the next in the path to v

        #region ComputeGeodesicDistance

        /// <summary>
        /// Calculate the geodesic distance, the shortest path
        /// on the surface, between this vertex and all other vertices
        /// uses Dijkstra's algorithm.
        /// The labels of vertices are set to point to the next in the path to v
        /// </summary>
        public void ComputeGeodesicDistance(Vertex v)
        {
            UnicityList<Vertex> set = new UnicityList<Vertex>();
            set.Add(v);
            Dijkstra(set);
        }

        public void ComputeGeodesicDistance(UnicityList<Vertex> set)
        {
            Dijkstra(set);
        }

        #endregion

        #region Dijkstra

        public void Dijkstra(UnicityList<Vertex> known)
        {
            int i, j;
            UnicityList<Vertex> S = new UnicityList<Vertex>();
            UnicityList<Vertex> W = new UnicityList<Vertex>();
            Vertex x, v;
            double min, d;

            for (j = 0; j < Vertices.Count; j++)
            {
                Vertices[j].label = -1;
                Vertices[j].Distance = 1e15;
                Vertices[j].Next = null;
            }
            //apply Dijkstra's algorithm
            for (j = 0; j < known.Count; j++)
            {
                v = known[j];
                v.Distance = 0;
                v.label = 0;
            }
            for (j = 0; j < known.Count; j++)
            {
                v = known[j];
                //set distance for neighbouring edges
                for (i = 0; i < v.Edges.Count; i++)
                {
                    d = v.Edges[i].Length;
                    if (v.Edges[i].Vertices[0] == v)
                    {
                        if (v.Edges[i].Vertices[1].Distance > d)
                        {
                            v.Edges[i].Vertices[1].Distance = d;
                            v.Edges[i].Vertices[1].Next = v;
                        }
                    }
                    else
                    {
                        if (v.Edges[i].Vertices[0].Distance > d)
                        {
                            v.Edges[i].Vertices[0].Distance = d;
                            v.Edges[i].Vertices[0].Next = v;
                        }
                    }
                }

                S.Add(v);

                //compute W, the set of all vertices connected to S, but not in S
                for (i = 0; i < v.Edges.Count; i++)
                {
                    if (v.Edges[i].Vertices[0] == v)
                    {
                        if (!S.Contains(v.Edges[i].Vertices[1]))
                            W.Add(v.Edges[i].Vertices[1]);
                    }
                    else
                    {
                        if (!S.Contains(v.Edges[i].Vertices[0]))
                            W.Add(v.Edges[i].Vertices[0]);
                    }
                }
            }

            while ((S.Count != Vertices.Count) && (W.Count > 0))
            {
                //find the vertex X in W for which the distance is minimum
                //there will be at least one element in W, because we stop when
                //   W is empty (and S contains all vertices);
                x = W[0];
                min = x.Distance;
                for (i = 1; i < W.Count; i++)
                {
                    if (W[i].Distance < min)
                    {
                        x = W[i];
                        min = x.Distance;
                    }
                }
                //Add X to S
                S.Add(x);
                //update W, the set of all vertices connected to S, but not in S
                W.Remove(x);
                //update closest distance for all vertices outside S (but only if they are connected to X), and update W
                for (i = 0; i < x.Edges.Count; i++)
                {
                    if (!(S.Contains(x.Edges[i].Vertices[0])))
                    {
                        d = x.Edges[i].Length;
                        if (d + x.Distance < x.Edges[i].Vertices[0].Distance)
                        {
                            x.Edges[i].Vertices[0].Distance = d + x.Distance;
                            //label prior vertex in sequence
                            x.Edges[i].Vertices[0].Next = x;
                        }
                        W.Add(x.Edges[i].Vertices[0]);
                    }
                    if (!(S.Contains(x.Edges[i].Vertices[1])))
                    {
                        d = x.Edges[i].Length;
                        if (d + x.Distance < x.Edges[i].Vertices[1].Distance)
                        {
                            x.Edges[i].Vertices[1].Distance = d + x.Distance;
                            //label prior vertex in sequence
                            x.Edges[i].Vertices[1].Next = x;
                        }
                        W.Add(x.Edges[i].Vertices[1]);
                    }
                }
                //next iteration
            }
        }

        #endregion

        public void DiscDeleteEdge(Edge e)
        {
            if (test(e.label & BOUNDARY)) return;
            e.label |= DELETED;
        }

        public void DiscDeleteFace(Face f, UnicityList<Edge> candidates)
        {
            int i;

            if (test(f.label & DELETED)) return;
            f.label |= DELETED;
            for (i = 0; i < 3; i++)
            {
                if (f.Edges[i].Faces[0] == f)
                {
                    f.Edges[i].label |= F0DEL;
                }
                else
                {
                    f.Edges[i].label |= F1DEL;
                }
                if ((test(f.Edges[i].label & F0DEL)) && (test(f.Edges[i].label & F1DEL)))
                {
                    //DiscDeleteEdge(f.edges[i]);
                }
                else
                {
                    candidates.Add(f.Edges[i]);
                }
            }
        }

        public Edge AddEdge(int v1, int v2)
        {
            int i;
            Vertex v;
            Edge e = new Edge();
            e.Vertices[0] = Vertices[v1];
            e.Vertices[1] = Vertices[v2];
            //if the edge exists, both vertices will have the edge
            //   in their adjacency lists
            v = Vertices[v1];
            for (i = 0; i < v.Edges.Count; i++)
            {
                if (e.Equals(v.Edges[i]))
                {
                    return v.Edges[i];
                }
            }
            Edges.AddWithoutCheck(e);
            Vertices[v1].Edges.Add(e);
            Vertices[v2].Edges.Add(e);
            return e;
        }

        public UnicityList<Vertex> ComputeRing(Vertex v, int level)
        {
            UnicityList<Vertex> ring = new UnicityList<Vertex>();
            ring.Add(v);
            if (level <= 0) return ring;
            return ComputeRing(ring, level);
        }

        public UnicityList<Vertex> ComputeRing(UnicityList<Vertex> set, int level)
        {
            UnicityList<Vertex> ring = new UnicityList<Vertex>();
            Vertex v;
            int i, j;

            for (i = 0; i < set.Count; i++)
            {
                v = set[i];
                ring.Add(v);
            }

            if (level <= 0) return ring;

            for (i = 0; i < set.Count; i++)
            {
                v = set[i];
                for (j = 0; j < v.Edges.Count; j++)
                {
                    ring.Add(v.Edges[j].Vertices[1]);
                    ring.Add(v.Edges[j].Vertices[0]);
                }
            }
            if (level == 1) return ring;
            return ComputeRing(ring, level - 1);
        }

        public UnicityList<Vertex> ComputeExclusiveRing(Vertex v, int level)
        {
            UnicityList<Vertex> ring = new UnicityList<Vertex>();
            if (level <= 0) return ring;
            ring.Add(v);
            return ComputeExclusiveRing(ring, level);
        }

        public UnicityList<Vertex> ComputeExclusiveRing(Face t, int level)
        {
            UnicityList<Vertex> ring = new UnicityList<Vertex>();
            if (level <= 0) return ring;
            ring.Add(t.Vertices[0]);
            ring.Add(t.Vertices[1]);
            ring.Add(t.Vertices[2]);
            return ComputeExclusiveRing(ring, level);
        }

        public UnicityList<Vertex> ComputeExclusiveRing(UnicityList<Vertex> set, int level)
        {
            UnicityList<Vertex> ring, inner;
            Vertex v;
            int i;

            inner = ComputeRing(set, level - 1);
            ring = ComputeRing(set, level);
            for (i = 0; i < inner.Count; i++)
            {
                v = inner[i];
                ring.Remove(v);
            }
            return ring;
        }

        public void CircleBoundary(UnicityList<Vertex> vboundary)
        {
            int i, n;
            Vector3 p = new Vector3();
            //Circle boundary
            n = vboundary.Count;
            for (i = 0; i < n; i++)
            {
                p.z = 0;
                p.x = System.Math.Cos(i * 2.0 * System.Math.PI / n) * 0.5 + 0.5;
                p.y = System.Math.Sin(i * 2.0 * System.Math.PI / n) * 0.5 + 0.5;
                vboundary[i].TexCoordinate = (Vector3)p.Clone();
            }
        }

        public void CircleChordBoundary(UnicityList<Vertex> vboundary)
        {
            int i, n;
            double dist, cd;
            Vector3 p = new Vector3();
            //Circle boundary
            n = vboundary.Count;
            dist = 0;
            for (i = 1; i < n; i++)
            {
                dist += (vboundary[i].Coordinate - vboundary[i - 1].Coordinate).Length;
            }
            dist += (vboundary[n - 1].Coordinate - vboundary[0].Coordinate).Length;
            cd = 0;
            for (i = 0; i < n; i++)
            {
                p.z = 0;
                p.x = System.Math.Cos(cd / dist * 2.0 * System.Math.PI) * 0.5 + 0.5;
                p.y = System.Math.Sin(cd / dist * 2.0 * System.Math.PI) * 0.5 + 0.5;
                vboundary[i].TexCoordinate = (Vector3)p.Clone();
                if (i < n - 1)
                    cd += (vboundary[i + 1].Coordinate - vboundary[i].Coordinate).Length;
            }
        }

        public void RectangleBoundary(UnicityList<Vertex> vboundary)
        {
            int i, n, basei, corner;
            double min;
            Vector3 p = new Vector3();
            //parameterise in a straight line
            n = vboundary.Count;
            if (n == 0) return;
            p.x = p.y = p.z = 0.0;
            for (i = 0; i < n; i++)
            {
                p.x = i / ((double)n);
                vboundary[i].TexCoordinate = (Vector3)p.Clone();
            }
            //find point closest to 1/4 distance, make it the corner
            basei = 0;
            corner = 0;
            min = System.Math.Abs(vboundary[basei].TexCoordinate.x - 1.0 / 4.0);
            for (i = basei; i < n; i++)
            {
                if (System.Math.Abs(vboundary[i].TexCoordinate.x - 1.0 / 4.0) < min)
                {
                    min = System.Math.Abs(vboundary[i].TexCoordinate.x - 1.0 / 4.0);
                    corner = i;
                }
            }
            for (i = basei; i < corner; i++)
            {
                p.x = p.z = 0.0;
                p.y = i / ((double)corner);
                vboundary[i].TexCoordinate = (Vector3)p.Clone();
            }

            //find point closest to 2/4 distance, make it the corner
            basei = corner;
            corner = basei;
            min = System.Math.Abs(vboundary[basei].TexCoordinate.x - 2.0 / 4.0);
            for (i = basei; i < n; i++)
            {
                if (System.Math.Abs(vboundary[i].TexCoordinate.x - 2.0 / 4.0) < min)
                {
                    min = System.Math.Abs(vboundary[i].TexCoordinate.x - 2.0 / 4.0);
                    corner = i;
                }
            }
            for (i = basei; i < corner; i++)
            {
                p.z = 0.0;
                p.y = 1.0;
                p.x = (i - basei) / ((double)(corner - basei));
                vboundary[i].TexCoordinate = (Vector3)p.Clone();
            }

            //find point closest to 3/4 distance, make it the corner
            basei = corner;
            corner = basei;
            min = System.Math.Abs(vboundary[basei].TexCoordinate.x - 3.0 / 4.0);
            for (i = basei; i < n; i++)
            {
                if (System.Math.Abs(vboundary[i].TexCoordinate.x - 3.0 / 4.0) < min)
                {
                    min = System.Math.Abs(vboundary[i].TexCoordinate.x - 3.0 / 4.0);
                    corner = i;
                }
            }
            for (i = basei; i < corner; i++)
            {
                p.z = 0.0;
                p.x = 1.0;
                p.y = 1.0 - (i - basei) / ((double)(corner - basei));
                vboundary[i].TexCoordinate = (Vector3)p.Clone();
            }

            //parameterise last stretch of boundary
            basei = corner;
            for (i = basei; i < n; i++)
            {
                p.z = 0.0;
                p.y = 0.0;
                p.x = 1.0 - (i - basei) / ((double)(n - basei));
                vboundary[i].TexCoordinate = (Vector3)p.Clone();
            }
        }

        public void RectangleChordBoundary(UnicityList<Vertex> vboundary)
        {
            int i, j, n, corner;
            double min;
            double dist, cd;
            Vector3 p = new Vector3();
            //parameterise in a straight line
            n = vboundary.Count;
            if (n == 0) return;
            p.x = p.y = p.z = 0.0;

            dist = 0;
            for (i = 1; i < n; i++)
            {
                dist += (vboundary[i].Coordinate - vboundary[i - 1].Coordinate).Length;
            }
            dist += (vboundary[n - 1].Coordinate - vboundary[0].Coordinate).Length;
            cd = 0;
            for (i = 0; i < n; i++)
            {
                p.x = cd / dist;
                vboundary[i].TexCoordinate = (Vector3)p.Clone();
                if (i < n - 1)
                    cd += (vboundary[i + 1].Coordinate - vboundary[i].Coordinate).Length;
            }

            //find corners
            //first corner is already at 0,0
            //last corner is the first corner
            for (j = 1; j < 4; j++)
            {
                //find closest parameter
                min = System.Math.Abs(vboundary[0].TexCoordinate.x - j / 4.0);
                corner = 0;
                for (i = 1; i < n; i++)
                {
                    if (System.Math.Abs(vboundary[i].TexCoordinate.x - j / 4.0) < min)
                    {
                        min = System.Math.Abs(vboundary[i].TexCoordinate.x - j / 4.0);
                        corner = i;
                    }
                }
                //set the parameter value
                vboundary[corner].TexCoordinate.x = j / 4.0;
            }
            //move each parameter to the correct segment
            for (i = 0; i < n; i++)
            {
                p = (Vector3)vboundary[i].TexCoordinate.Clone();
                //first segment
                if (p.x < 1.0 / 4.0)
                {
                    p.y = p.x / 0.25;
                    p.x = p.z = 0.0;
                }
                else if (p.x < 2.0 / 4.0)
                {
                    p.x = (p.x - 1.0 / 4.0) / 0.25;
                    p.y = 1.0;
                    p.z = 0.0;
                }
                else if (p.x < 3.0 / 4.0)
                {
                    p.y = 1.0 - (p.x - 2.0 / 4.0) / 0.25;
                    p.x = 1.0;
                    p.z = 0.0;
                }
                else
                {
                    p.x = 1.0 - (p.x - 3.0 / 4.0) / 0.25;
                    p.y = 0.0;
                    p.z = 0.0;
                }
                vboundary[i].TexCoordinate = (Vector3)p.Clone();
            }
        }

        public void ComputeBoundary(UnicityList<Edge> eboundary)
        {
            int i;

            //mark boundary edges
            for (i = 0; i < Edges.Count; i++)
            {
                if (Edges[i].Faces.Count < 2)
                {
                    eboundary.Add(Edges[i]);
                }
            }
        }

        //Only works if there is precisely one boundary loop!!
        public void ComputeBoundary(UnicityList<Edge> eboundary, UnicityList<Vertex> vboundary)
        {
            int i;
            Edge e;
            Vertex v;

            //mark boundary edges
            Console.Write("Total edges: {0}\n", Edges.Count);
            for (i = 0; i < Edges.Count; i++)
            {
                Edges[i].label = 0;
                //mark boundaries
                if (Edges[i].Faces.Count < 2)
                {
                    Edges[i].label |= BOUNDARY;
                    eboundary.Add(Edges[i]);
                }
            }
            Console.Write("Boundary edges: {0}\n", eboundary.Count);

            //remove edges one for one to get boundary vertices
            //srand(time(null));
            //e=eboundary[rand()%eboundary.getCount()];
            if (eboundary.Count == 0) Console.Write("Fatal: bounray count=0\n");
            e = eboundary[0];
            v = e.Vertices[1];
            vboundary.Add(e.Vertices[0]);
            vboundary.Add(e.Vertices[1]);
            eboundary.Remove(0);

            while (eboundary.Count > 0)
            {
                for (i = 0; i < eboundary.Count; i++)
                {
                    e = eboundary[i];
                    if (e.Vertices[0] == v)
                    {
                        vboundary.Add(e.Vertices[1]);
                        eboundary.Remove(i);
                        v = e.Vertices[1];
                        e = null;
                        break;
                    }
                    else if (e.Vertices[1] == v)
                    {
                        vboundary.Add(e.Vertices[0]);
                        eboundary.Remove(i);
                        v = e.Vertices[0];
                        e = null;
                        break;
                    }
                }
                if (e != null)
                    Console.Write("Broken loop! {0} edges left\n", eboundary.Count);
                if (e != null)
                    break;
            }

        }

        //Tutte
        public void ParamBarycentric(int boundary_type)
        {
            int i, j, n, c;
            double lambda;
            UnicityList<Edge> eboundary = new UnicityList<Edge>();
            UnicityList<Vertex> vboundary = new UnicityList<Vertex>();
            Vertex v, v2;
            SparseMatrix M = new SparseMatrix();
            double[] bu;
            double[] bv;
            //for second chance to solve
            double[] scu;
            double[] scv;

            ComputeBoundary(eboundary, vboundary);
            if (vboundary.Count == 0)
            {
                Console.Write("No boundary, perhaps the mesh is too small?\n");
                return;
            }
            switch (boundary_type)
            {
                case BOUNDARY_CIRCLE:
                    CircleBoundary(vboundary);
                    break;
                case BOUNDARY_RECTANGLE:
                    RectangleBoundary(vboundary);
                    break;
                case BOUNDARY_CHORD_CIRCLE:
                    CircleChordBoundary(vboundary);
                    break;
                case BOUNDARY_CHORD_RECTANGLE:
                    RectangleChordBoundary(vboundary);
                    break;
                case BOUNDARY_PRESET:
                    break;
                default:
                    RectangleBoundary(vboundary);
                    break;
            }

            //label non-border vertices according to equation variable number
            n = 0;
            for (i = 0; i < Vertices.Count; i++)
            {
                if (!vboundary.Contains(Vertices[i]))
                {
                    Vertices[i].label = n;
                    n++;
                }
                else Vertices[i].label = -1;
            }

            M.Resize(n, n);
            bu = new double[n];
            for (i = 0; i < n; i++) bu[i] = 0.0;
            bv = new double[n];
            for (i = 0; i < n; i++) bv[i] = 0.0;
            scu = new double[n];
            scv = new double[n];
            //border is paramerised. Set up system of equations setting each vertex to the
            //   barycentre of its neighbours (Tutte)
            for (i = 0; i < Vertices.Count; i++)
            {
                //check for a boundary
                v = Vertices[i];
                if (v.label != -1)
                {
                    c = v.label;
                    M[c, c] = 1.0;
                    //lambda=1.0/v.valence();
                    //Compute the 1 ring
                    UnicityList<Vertex> ring;
                    ring = ComputeExclusiveRing(v, 1);
                    //From Floater: lambda=w_ij / sum(w_ij)
                    // for barycentric w_ij=1
                    lambda = 1.0 / (ring.Count);

                    for (j = 0; j < ring.Count; j++)
                    {
                        v2 = ring[j];
                        if (v2 == v) continue;
                        //boundary goes into b
                        if (v2.label == -1)
                        {
                            bu[c] += lambda * v2.TexCoordinate.x;
                            bv[c] += lambda * v2.TexCoordinate.y;
                        }
                        else
                        {
                            M[c, v2.label] = -lambda;
                        }
                    }
                }
            }

            for (i = 0; i < n; i++) scu[i] = bu[i];
            for (i = 0; i < n; i++) scv[i] = bv[i];
            c = M.GaussSeidel(bu, 100000, 1e-7);
            if (c < 0)
            {
                Console.Write("Failed to solve\n");
                return;
            }
            if (c == 0)
            {
                //failed to obtain accuracy, try again
                c = M.GaussSeidel(scu, 100000, 1e-7, bu);
                if (c == 0)
                    Console.Write("Failed to solve to desired accuracy\n");
                for (i = 0; i < n; i++) bu[i] = scu[i];
            }
            c = M.GaussSeidel(bv, 100000, 1e-7);
            if (c < 0)
            {
                Console.Write("Failed to solve\n");
                return;
            }
            if (c == 0)
            {
                //failed to obtain accuracy, try again
                c = M.GaussSeidel(scv, 100000, 1e-7, bv);
                if (c == 0)
                    Console.Write("Failed to solve to desired accuracy\n");
                for (i = 0; i < n; i++) bv[i] = scv[i];
            }

            for (i = 0; i < Vertices.Count; i++)
            {
                //check for a boundary
                v = Vertices[i];
                if (v.label != -1)
                {
                    c = v.label;
                    v.TexCoordinate.x = bu[c];
                    v.TexCoordinate.y = bv[c];
                    v.TexCoordinate.z = 0.0;
                }
            }
        }


        //Stretch for a vertex
        public double Sigma(Vertex v)
        {
            int i;
            double sumarea = 0.0, sumstretch = 0.0;
            double area, stretch;

            for (i = 0; i < v.Faces.Count; i++)
            {
                area = v.Faces[i].Area();
                sumarea += area;
                stretch = v.Faces[i].Stretch();
                sumstretch += area * stretch * stretch;
            }

            return System.Math.Sqrt(sumstretch / sumarea);
        }

        //Stretch for entire mesh
        public double AllStretch()
        {
            int i;
            double sumarea = 0.0, sumstretch = 0.0;
            double area, stretch;

            for (i = 0; i < Faces.Count; i++)
            {
                area = Faces[i].Area();
                sumarea += area;
                stretch = Faces[i].Stretch();
                sumstretch += area * stretch * stretch;
            }

            return System.Math.Sqrt(sumstretch / sumarea);
        }

        //area of triangle
        public double area(Vector3 v1, Vector3 v2, Vector3 v3)
        {
            return (v2 - v1).Cross(v3 - v1).Length / 2.0;
        }

        public void ParamFloater(int boundary_type)
        {
            int i, j, k, n, c;
            double lambda, sumtheta, theta, totalw;
            UnicityList<Edge> eboundary = new UnicityList<Edge>();
            UnicityList<Vertex> vboundary = new UnicityList<Vertex>();
            Vertex v, v2 = null, v3 = null;
            SparseMatrix M = new SparseMatrix();
            double[] bu;
            double[] bv;
            //for second chance to solve
            double[] scu;
            double[] scv;
            bool trianglesleft = true;

            Console.Write("Boundary\n");
            ComputeBoundary(eboundary, vboundary);
            if (vboundary.Count == 0)
            {
                Console.Write("No boundary, perhaps the mesh is too small?\n");
                return;
            }
            Console.Write("Boundary param\n");
            switch (boundary_type)
            {
                case BOUNDARY_CIRCLE:
                    CircleBoundary(vboundary);
                    break;
                case BOUNDARY_RECTANGLE:
                    RectangleBoundary(vboundary);
                    break;
                case BOUNDARY_CHORD_CIRCLE:
                    CircleChordBoundary(vboundary);
                    break;
                case BOUNDARY_CHORD_RECTANGLE:
                    RectangleChordBoundary(vboundary);
                    break;
                case BOUNDARY_PRESET:
                    break;
                default:
                    RectangleBoundary(vboundary);
                    break;
            }

            //label non-border vertices according to equation variable number
            n = 0;
            for (i = 0; i < Vertices.Count; i++)
            {
                if (!vboundary.Contains(Vertices[i]))
                {
                    Vertices[i].label = n;
                    n++;
                }
                else Vertices[i].label = -1;
            }

            Console.Write("Main param...\n");
            M.Resize(n, n);
            bu = new double[n];
            for (i = 0; i < n; i++) bu[i] = 0.0;
            bv = new double[n];
            for (i = 0; i < n; i++) bv[i] = 0.0;
            scu = new double[n];
            scv = new double[n];
            //border is paramerised. Set up system of equations setting each vertex to the
            //   barycentre of its neighbours (Tutte)
            for (i = 0; i < Vertices.Count; i++)
            {
                //check for a boundary
                v = Vertices[i];
                //corresponds to edges
                v.wij = new double[v.Edges.Count];
                for (j = 0; j < v.Edges.Count; j++)
                {
                    v.wij[j] = 0.0;
                }

                if (v.label != -1)
                {
                    UnicityList<Vertex> ring;
                    ring = ComputeExclusiveRing(v, 1);
                    //save texture coordinates of boundary in normals
                    for (j = 0; j < ring.Count; j++)
                    {
                        ring[j].n = (Vector3)ring[j].TexCoordinate.Clone();
                    }
                    c = v.label;
                    M[c, c] = 1.0;
                    //Now do local parameterisation according to Floater
                    //set p in the middle
                    v.TexCoordinate = new Vector3(0, 0, 0);
                    //determine sum of angles 
                    sumtheta = 0.0;
                    for (j = 0; j < v.Faces.Count; j++)
                    {
                        if (v == v.Faces[j].Vertices[0])
                        {
                            v2 = v.Faces[j].Vertices[1];
                            v3 = v.Faces[j].Vertices[2];
                        }
                        else if (v == v.Faces[j].Vertices[1])
                        {
                            v2 = v.Faces[j].Vertices[0];
                            v3 = v.Faces[j].Vertices[2];
                        }
                        else if (v == v.Faces[j].Vertices[2])
                        {
                            v2 = v.Faces[j].Vertices[0];
                            v3 = v.Faces[j].Vertices[1];
                        }
                        sumtheta += System.Math.Abs((v2.Coordinate - v.Coordinate).Angle(v3.Coordinate - v.Coordinate));
                        //this face has not been parameterised yet
                        v.Faces[j].label = 0;
                    }

                    //place first adjacent vertex horizontally, preserving distance
                    j = 0;
                    if (v == v.Faces[j].Vertices[0])
                    {
                        v2 = v.Faces[j].Vertices[1];
                        v3 = v.Faces[j].Vertices[2];
                    }
                    else if (v == v.Faces[j].Vertices[1])
                    {
                        v2 = v.Faces[j].Vertices[0];
                        v3 = v.Faces[j].Vertices[2];
                    }
                    else if (v == v.Faces[j].Vertices[2])
                    {
                        v2 = v.Faces[j].Vertices[0];
                        v3 = v.Faces[j].Vertices[1];
                    }
                    theta = 0.0;
                    theta += 2 * System.Math.PI * System.Math.Abs((v2.Coordinate - v.Coordinate).Angle(v3.Coordinate - v.Coordinate)) / sumtheta;
                    v2.TexCoordinate.x = (v2.Coordinate - v.Coordinate).Length;
                    v2.TexCoordinate.y = v2.TexCoordinate.z = 0.0;
                    v3.TexCoordinate.x = System.Math.Cos(theta);
                    v3.TexCoordinate.y = System.Math.Sin(theta);
                    v3.TexCoordinate.z = 0.0;
                    v3.TexCoordinate *= (v3.Coordinate - v.Coordinate).Length;
                    v2 = v3;

                    //place vertices at shape preserving angles around the vertex
                    //this creates a local shape preserving parameterization
                    trianglesleft = true;
                    while (trianglesleft)
                    {
                        //find a triangle adjacent to the current one, that 
                        //has not been computed 
                        Edge e = new Edge();
                        e.Vertices[0] = v;
                        e.Vertices[1] = v2;

                        trianglesleft = false;
                        for (j = 0; j < v.Faces.Count; j++)
                        {
                            if (v.Faces[j].label == 0)
                            {
                                //check to see if common edge is shared
                                if (v.Faces[j].Edges[0].Equals(e)) break;
                                if (v.Faces[j].Edges[1].Equals(e)) break;
                                if (v.Faces[j].Edges[2].Equals(e)) break;
                            }
                        }
                        //found an adjacent triangle
                        if (j < v.Faces.Count)
                        {
                            int[] tag = new int[3] { 0, 0, 0 };
                            for (k = 0; k < 3; k++)
                            {
                                if (v == v.Faces[j].Vertices[k]) tag[k] = 1;
                                if (v2 == v.Faces[j].Vertices[k]) tag[k] = 1;
                            }
                            for (k = 0; k < 3; k++)
                            {
                                if (tag[k] == 0) v3 = v.Faces[j].Vertices[k];
                            }
                            //the undetermined vertex is placed according to
                            //its angle to the vertex, and the parametric distance
                            //   is set to be the same as the real distance
                            theta += 2 * System.Math.PI * System.Math.Abs((v2.Coordinate - v.Coordinate).Angle(v3.Coordinate - v.Coordinate)) / sumtheta;
                            v3.TexCoordinate.x = System.Math.Cos(theta);
                            v3.TexCoordinate.y = System.Math.Sin(theta);
                            v3.TexCoordinate.z = 0.0;
                            v3.TexCoordinate *= (v3.Coordinate - v.Coordinate).Length;
                            v2 = v3;
                            trianglesleft = true;
                            v.Faces[j].label = 1;
                        }
                    }

                    //for each vertex determine a triangle that  
                    //  v lies in, and then adjust the weights
                    //  according to the barycentric coordinates
                    //Note: ring.getCount()==v.valence()
                    for (k = 0; k < ring.Count; k++)
                    {
                        for (j = 0; j < v.Faces.Count; j++)
                        {
                            if (v == v.Faces[j].Vertices[0])
                            {
                                v2 = v.Faces[j].Vertices[1];
                                v3 = v.Faces[j].Vertices[2];
                            }
                            else if (v == v.Faces[j].Vertices[1])
                            {
                                v2 = v.Faces[j].Vertices[0];
                                v3 = v.Faces[j].Vertices[2];
                            }
                            else if (v == v.Faces[j].Vertices[2])
                            {
                                v2 = v.Faces[j].Vertices[0];
                                v3 = v.Faces[j].Vertices[1];
                            }
                            //determine if this is the triangle
                            if (Vector3.IsCrossing(ring[k].TexCoordinate, v.TexCoordinate, v2.TexCoordinate, v3.TexCoordinate))
                            {
                                double l1, l2, l3, a;
                                int i2, i3;
                                i2 = ring.IndexOf(v2);
                                i3 = ring.IndexOf(v3);
                                if (i2 < 0) Console.Write("Bad index: should not happen!\n");
                                if (i3 < 0) Console.Write("Bad index: should not happen!\n");
                                a = area(ring[k].TexCoordinate, v2.TexCoordinate, v3.TexCoordinate);
                                l1 = area(v.TexCoordinate, v2.TexCoordinate, v3.TexCoordinate) / a;
                                l2 = area(ring[k].TexCoordinate, v.TexCoordinate, v3.TexCoordinate) / a;
                                l3 = area(ring[k].TexCoordinate, v2.TexCoordinate, v.TexCoordinate) / a;
                                v.wij[k] += l1;
                                v.wij[i2] += l2;
                                v.wij[i3] += l3;
                            }
                        }
                    }

                    totalw = 0.0;
                    for (j = 0; j < ring.Count; j++)
                    {
                        totalw += v.wij[j];
                        //if it is a border vertex, put the normal
                        // back in the texture
                        if (ring[j].label == -1)
                        {
                            ring[j].TexCoordinate = (Vector3)ring[j].n.Clone();
                        }
                    }
                    //if (totalw==0.0) totalw=1.0;

                    for (j = 0; j < ring.Count; j++)
                    {
                        v2 = ring[j];
                        lambda = v.wij[j] / totalw;
                        if (v2 == v) continue;
                        //boundary goes into b
                        if (v2.label == -1)
                        {
                            bu[c] += lambda * v2.TexCoordinate.x;
                            bv[c] += lambda * v2.TexCoordinate.y;
                        }
                        else
                        {
                            M[c, v2.label] = -lambda;
                        }
                    }
                }
            }

            for (i = 0; i < n; i++) scu[i] = bu[i];
            for (i = 0; i < n; i++) scv[i] = bv[i];
            Console.Write("Solving\n");
            c = M.GaussSeidel(bu, 100000, 1e-7);
            //c=M.ConjugateGradient(bu, 100000, 1e-7);
            if (c < 0)
            {
                Console.Write("Failed to solve\n");
                return;
            }
            if (c == 0)
            {
                //failed to obtain accuracy, try again
                c = M.GaussSeidel(scu, 100000, 1e-7, bu);
                //c=M.ConjugateGradient(scu, 100000, 1e-7, bu);
                if (c == 0)
                    Console.Write("Failed to solve to desired accuracy\n");
                for (i = 0; i < n; i++) bu[i] = scu[i];
            }
            c = M.GaussSeidel(bv, 100000, 1e-7);
            //c=M.ConjugateGradient(bv, 100000, 1e-7);
            if (c < 0)
            {
                Console.Write("Failed to solve\n");
                return;
            }
            if (c == 0)
            {
                //failed to obtain accuracy, try again
                c = M.GaussSeidel(scv, 100000, 1e-7, bv);
                //c=M.ConjugateGradient(scv, 100000, 1e-7, bv);
                if (c == 0)
                    Console.Write("Failed to solve to desired accuracy\n");
                for (i = 0; i < n; i++) bv[i] = scv[i];
            }

            for (i = 0; i < Vertices.Count; i++)
            {
                //check for a boundary
                v = Vertices[i];
                if (v.label != -1)
                {
                    c = v.label;
                    v.TexCoordinate.x = bu[c];
                    v.TexCoordinate.y = bv[c];
                    v.TexCoordinate.z = 0.0;
                }
            }
        }

        //Mean value parameterisation, also due to Floater
        public void ParamMeanValue(int boundary_type)
        {
            int i, j, k, n, c;
            double lambda, sumtheta, theta, totalw;
            UnicityList<Edge> eboundary = new UnicityList<Edge>();
            UnicityList<Vertex> vboundary = new UnicityList<Vertex>();
            Vertex v, v2 = null, v3 = null;
            SparseMatrix M = new SparseMatrix();
            double[] bu;
            double[] bv;
            //for second chance to solve
            double[] scu;
            double[] scv;
            bool trianglesleft = true;

            ComputeBoundary(eboundary, vboundary);
            if (vboundary.Count == 0)
            {
                Console.Write("No boundary, perhaps the mesh is too small?\n");
                return;
            }
            switch (boundary_type)
            {
                case BOUNDARY_CIRCLE:
                    CircleBoundary(vboundary);
                    break;
                case BOUNDARY_RECTANGLE:
                    RectangleBoundary(vboundary);
                    break;
                case BOUNDARY_CHORD_CIRCLE:
                    CircleChordBoundary(vboundary);
                    break;
                case BOUNDARY_CHORD_RECTANGLE:
                    RectangleChordBoundary(vboundary);
                    break;
                case BOUNDARY_PRESET:
                    break;
                default:
                    RectangleBoundary(vboundary);
                    break;
            }

            //label non-border vertices according to equation variable number
            n = 0;
            for (i = 0; i < Vertices.Count; i++)
            {
                if (!vboundary.Contains(Vertices[i]))
                {
                    Vertices[i].label = n;
                    n++;
                }
                else Vertices[i].label = -1;
            }

            M.Resize(n, n);
            bu = new double[n];
            for (i = 0; i < n; i++) bu[i] = 0.0;
            bv = new double[n];
            for (i = 0; i < n; i++) bv[i] = 0.0;
            scu = new double[n];
            scv = new double[n];
            //border is paramerised. Set up system of equations setting each vertex to the
            //   barycentre of its neighbours (Tutte)
            for (i = 0; i < Vertices.Count; i++)
            {
                //check for a boundary
                v = Vertices[i];
                //corresponds to edges
                v.wij = new double[v.Edges.Count];
                for (j = 0; j < v.Edges.Count; j++)
                {
                    v.wij[j] = 0.0;
                }

                if (v.label != -1)
                {
                    UnicityList<Vertex> ring;
                    ring = ComputeExclusiveRing(v, 1);
                    //save texture coordinates of boundary in normals
                    for (j = 0; j < ring.Count; j++)
                    {
                        ring[j].Normal = (Vector3)ring[j].TexCoordinate.Clone();
                    }
                    c = v.label;
                    M[c, c] = 1.0;
                    //Now do local parameterisation according to Floater
                    //set p in the middle
                    v.TexCoordinate = new Vector3(0, 0, 0);
                    //determine sum of angles 
                    sumtheta = 0.0;
                    for (j = 0; j < v.Faces.Count; j++)
                    {
                        if (v == v.Faces[j].Vertices[0])
                        {
                            v2 = v.Faces[j].Vertices[1];
                            v3 = v.Faces[j].Vertices[2];
                        }
                        else if (v == v.Faces[j].Vertices[1])
                        {
                            v2 = v.Faces[j].Vertices[0];
                            v3 = v.Faces[j].Vertices[2];
                        }
                        else if (v == v.Faces[j].Vertices[2])
                        {
                            v2 = v.Faces[j].Vertices[0];
                            v3 = v.Faces[j].Vertices[1];
                        }
                        sumtheta += System.Math.Abs((v2.Coordinate - v.Coordinate).Angle(v3.Coordinate - v.Coordinate));
                        //this face has not been parameterised yet
                        v.Faces[j].label = 0;
                    }

                    //place first adjacent vertex horizontally, preserving distance
                    j = 0;
                    if (v == v.Faces[j].Vertices[0])
                    {
                        v2 = v.Faces[j].Vertices[1];
                        v3 = v.Faces[j].Vertices[2];
                    }
                    else if (v == v.Faces[j].Vertices[1])
                    {
                        v2 = v.Faces[j].Vertices[0];
                        v3 = v.Faces[j].Vertices[2];
                    }
                    else if (v == v.Faces[j].Vertices[2])
                    {
                        v2 = v.Faces[j].Vertices[0];
                        v3 = v.Faces[j].Vertices[1];
                    }
                    theta = 0.0;
                    theta += 2 * System.Math.PI * System.Math.Abs((v2.Coordinate - v.Coordinate).Angle(v3.Coordinate - v.Coordinate)) / sumtheta;
                    v2.TexCoordinate.x = (v2.Coordinate - v.Coordinate).Length;
                    v2.TexCoordinate.y = v2.TexCoordinate.z = 0.0;
                    v3.TexCoordinate.x = System.Math.Cos(theta);
                    v3.TexCoordinate.y = System.Math.Sin(theta);
                    v3.TexCoordinate.z = 0.0;
                    v3.TexCoordinate *= (v3.Coordinate - v.Coordinate).Length;
                    v2 = v3;

                    //place vertices at shape preserving angles around the vertex
                    //this creates a local shape preserving parameterization
                    trianglesleft = true;
                    while (trianglesleft)
                    {
                        //find a triangle adjacent to the current one, that 
                        //has not been computed 
                        Edge e = new Edge();
                        e.Vertices[0] = v;
                        e.Vertices[1] = v2;

                        trianglesleft = false;
                        for (j = 0; j < v.Faces.Count; j++)
                        {
                            if (v.Faces[j].label == 0)
                            {
                                //check to see if common edge is shared
                                if (v.Faces[j].Edges[0].Equals(e)) break;
                                if (v.Faces[j].Edges[1].Equals(e)) break;
                                if (v.Faces[j].Edges[2].Equals(e)) break;
                            }
                        }
                        //found an adjacent triangle
                        if (j < v.Faces.Count)
                        {
                            int[] tag = new int[3] { 0, 0, 0 };
                            for (k = 0; k < 3; k++)
                            {
                                if (v == v.Faces[j].Vertices[k]) tag[k] = 1;
                                if (v2 == v.Faces[j].Vertices[k]) tag[k] = 1;
                            }
                            for (k = 0; k < 3; k++)
                            {
                                if (tag[k] == 0) v3 = v.Faces[j].Vertices[k];
                            }
                            //the undetermined vertex is placed according to
                            //its angle to the vertex, and the parametric distance
                            //   is set to be the same as the real distance
                            theta += 2 * System.Math.PI * System.Math.Abs((v2.Coordinate - v.Coordinate.Angle(v3.Coordinate - v.Coordinate)) / sumtheta;
                            v3.TexCoordinate.x = System.Math.Cos(theta);
                            v3.TexCoordinate.y = System.Math.Sin(theta);
                            v3.TexCoordinate.z = 0.0;
                            v3.TexCoordinate *= (v3.Coordinate - v.Coordinate).Length;
                            v2 = v3;
                            trianglesleft = true;
                            v.Faces[j].label = 1;
                        }
                    }

                    for (k = 0; k < ring.Count; k++)
                    {
                        double alpha, beta;
                        Edge e;
                        Edge e2 = new Edge();

                        e2.Vertices[0] = v;
                        e2.Vertices[1] = v2;
                        v2 = ring[k];
                        e = v.Edges.Match(e2);
                        //find angles adjacent to edge
                        Face f = e.Faces[0];
                        if (v == f.Vertices[0])
                        {
                            v2 = f.Vertices[1];
                            v3 = f.Vertices[2];
                        }
                        else if (v == f.Vertices[1])
                        {
                            v2 = f.Vertices[0];
                            v3 = f.Vertices[2];
                        }
                        else if (v == f.Vertices[2])
                        {
                            v2 = f.Vertices[0];
                            v3 = f.Vertices[1];
                        }
                        alpha = (v3.Coordinate - v.Coordinate).Angle(v2.Coordinate - v.Coordinate);

                        f = e.Faces[1];
                        if (v == f.Vertices[0])
                        {
                            v2 = f.Vertices[1];
                            v3 = f.Vertices[2];
                        }
                        else if (v == f.Vertices[1])
                        {
                            v2 = f.Vertices[0];
                            v3 = f.Vertices[2];
                        }
                        else if (v == f.Vertices[2])
                        {
                            v2 = f.Vertices[0];
                            v3 = f.Vertices[1];
                        }
                        beta = (v3.Coordinate - v.Coordinate.Angle(v2.Coordinate - v.Coordinate);

                        //Formula from mean-value paper by Floater
                        v.wij[k] = (System.Math.Tan(alpha / 2.0) + System.Math.Tan(beta / 2.0)) / (v2.Coordinate - v.Coordinate).Length;
                    }

                    totalw = 0.0;
                    for (j = 0; j < ring.Count; j++)
                    {
                        totalw += v.wij[j];
                        //if it is a border vertex, put the normal
                        // back in the texture
                        if (ring[j].label == -1)
                        {
                            ring[j].TexCoordinate = (Vector3)ring[j].Normal.Clone();
                        }
                    }
                    //if (totalw==0.0) totalw=1.0;

                    for (j = 0; j < ring.Count; j++)
                    {
                        v2 = ring[j];
                        lambda = v.wij[j] / totalw;
                        if (v2 == v) continue;
                        //boundary goes into b
                        if (v2.label == -1)
                        {
                            bu[c] += lambda * v2.TexCoordinate.x;
                            bv[c] += lambda * v2.TexCoordinate.y;
                        }
                        else
                        {
                            M[c, v2.label] = -lambda;
                        }
                    }
                }
            }

            for (i = 0; i < n; i++) scu[i] = bu[i];
            for (i = 0; i < n; i++) scv[i] = bv[i];
            c = M.GaussSeidel(bu, 100000, 1e-7);
            if (c < 0)
            {
                Console.Write("Failed to solve\n");
                return;
            }
            if (c == 0)
            {
                //failed to obtain accuracy, try again
                c = M.GaussSeidel(scu, 100000, 1e-7, bu);
                if (c == 0)
                    Console.Write("Failed to solve to desired accuracy\n");
                for (i = 0; i < n; i++) bu[i] = scu[i];
            }
            c = M.GaussSeidel(bv, 100000, 1e-7);
            if (c < 0)
            {
                Console.Write("Failed to solve\n");
                return;
            }
            if (c == 0)
            {
                //failed to obtain accuracy, try again
                c = M.GaussSeidel(scv, 100000, 1e-7, bv);
                if (c == 0)
                    Console.Write("Failed to solve to desired accuracy\n");
                for (i = 0; i < n; i++) bv[i] = scv[i];
            }

            for (i = 0; i < Vertices.Count; i++)
            {
                //check for a boundary
                v = Vertices[i];
                if (v.label != -1)
                {
                    c = v.label;
                    v.TexCoordinate.x = bu[c];
                    v.TexCoordinate.y = bv[c];
                    v.TexCoordinate.z = 0.0;
                }
            }
        }

        //Use Floater parameterisation, followed by stretch minimisation of Yoshizawa, Belyaev and Seidel (2004)
        public void ParamStretchmin(int boundary_type, double eta)
        {
            int i, j, n, c;
            double lambda, sigma, totalw, Enew, Eold;
            UnicityList<Edge> eboundary = new UnicityList<Edge>();
            UnicityList<Vertex> vboundary = new UnicityList<Vertex>();
            Vertex v, v2;
            SparseMatrix M = new SparseMatrix();
            double[] bu;
            double[] bv;
            double[] su;
            double[] sv;
            //for second chance to solve
            double[] scu;
            double[] scv;
            int iter = 500;
            Mesh opt = null;

            ParamMeanValue(boundary_type);
            //ParamBarycentric(boundary_type);
            ComputeBoundary(eboundary, vboundary);
            if (vboundary.Count == 0)
            {
                Console.Write("No boundary, perhaps the mesh is too small?\n");
                return;
            }
            switch (boundary_type)
            {
                case BOUNDARY_CIRCLE:
                    CircleBoundary(vboundary);
                    break;
                case BOUNDARY_RECTANGLE:
                    RectangleBoundary(vboundary);
                    break;
                case BOUNDARY_CHORD_CIRCLE:
                    CircleChordBoundary(vboundary);
                    break;
                case BOUNDARY_CHORD_RECTANGLE:
                    RectangleChordBoundary(vboundary);
                    break;
                case BOUNDARY_PRESET:
                    break;
                default:
                    RectangleBoundary(vboundary);
                    break;
            }

            //label non-border vertices according to equation variable number
            n = 0;
            for (i = 0; i < Vertices.Count; i++)
            {
                if (!vboundary.Contains(Vertices[i]))
                {
                    Vertices[i].label = n;
                    n++;
                }
                else Vertices[i].label = -1;
            }
            Console.Write("variables={0}\n", n);

            M.Resize(n, n);
            bu = new double[n];
            for (i = 0; i < n; i++) bu[i] = 0.0;
            bv = new double[n];
            for (i = 0; i < n; i++) bv[i] = 0.0;
            su = new double[n];
            sv = new double[n];
            scu = new double[n];
            scv = new double[n];

            Eold = AllStretch() + 1.0;
            while (--iter > 0)
            {
                Enew = AllStretch();
                Console.Write("Eold={0} Enew={1}\n", Eold, Enew);
                if ((Enew >= Eold) || (Double.IsNaN(Enew))) break;
                Eold = Enew;
                opt = this.Clone();
                Console.Write("Stretchmin iter {0}\n", iter);

                for (i = 0; i < n; i++) bu[i] = 0.0;
                for (i = 0; i < n; i++) bv[i] = 0.0;
                for (i = 0; i < Vertices.Count; i++)
                {
                    //check for a boundary
                    v = Vertices[i];

                    if (v.label != -1)
                    {
                        c = v.label;
                        M[c, c] = 1.0;
                        UnicityList<Vertex> ring;
                        ring = ComputeExclusiveRing(v, 1);
                        for (j = 0; j < ring.Count; j++)
                        {
                            sigma = Sigma(ring[j]);
                            sigma = System.Math.Pow(sigma, eta);
                            if (sigma != 0.0)
                                v.wij[j] /= sigma;
                            else
                                Console.Write("bad sigma\n");
                        }

                        totalw = 0.0;
                        for (j = 0; j < ring.Count; j++)
                        {
                            totalw += v.wij[j];
                        }
                        //if (totalw==0.0) {Console.Write("smallw\n"); totalw=1.0;}

                        for (j = 0; j < ring.Count; j++)
                        {
                            v2 = ring[j];
                            lambda = v.wij[j] / totalw;
                            if (lambda < 0) Console.Write("-lambda\n");
                            if (v2 == v) continue;
                            //boundary goes into b
                            if (v2.label == -1)
                            {
                                if (v2.TexCoordinate.x < 0.0) Console.Write("smallx\n");
                                if (v2.TexCoordinate.x > 1.0) Console.Write("largex\n");
                                if (v2.TexCoordinate.y < 0.0) Console.Write("smally\n");
                                if (v2.TexCoordinate.y > 1.0) Console.Write("largey\n");
                                bu[c] += lambda * v2.TexCoordinate.x;
                                bv[c] += lambda * v2.TexCoordinate.y;
                            }
                            else
                            {
                                M[c, v2.label] = -lambda;
                            }
                        }
                    }
                }

                for (i = 0; i < Vertices.Count; i++)
                {
                    //check for a boundary
                    v = Vertices[i];
                    if (v.label != -1)
                    {
                        c = v.label;
                        su[c] = v.TexCoordinate.x;
                        sv[c] = v.TexCoordinate.y;
                    }
                }

                for (i = 0; i < n; i++) scu[i] = bu[i];
                for (i = 0; i < n; i++) scv[i] = bv[i];
                //c=M.GaussSeidel(bu, 100000, 1e-7, su);
                c = M.GaussSeidel(bu, 100000, 1e-7);
                if (c < 0)
                {
                    Console.Write("Failed to solve\n");
                    Copy(opt);
                    return;
                }
                if (c == 0)
                {
                    //failed to obtain accuracy, try again
                    c = M.GaussSeidel(scu, 100000, 1e-7, bu);
                    if (c == 0)
                    { Console.Write("Failed to solve to desired accuracy\n"); eta /= 2.0; continue; }
                    for (i = 0; i < n; i++) bu[i] = scu[i];
                }
                //c=M.GaussSeidel(bv, 100000, 1e-7, sv);
                c = M.GaussSeidel(bv, 100000, 1e-7);
                if (c < 0)
                {
                    Console.Write("Failed to solve\n");
                    Copy(opt);
                    return;
                }
                if (c == 0)
                {
                    //failed to obtain accuracy, try again
                    c = M.GaussSeidel(scv, 100000, 1e-7, bv);
                    if (c == 0)
                    { Console.Write("Failed to solve to desired accuracy\n"); eta /= 2.0; continue; }
                    for (i = 0; i < n; i++) bv[i] = scv[i];
                }
                for (i = 0; i < Vertices.Count; i++)
                {
                    //check for a boundary
                    v = Vertices[i];
                    if (v.label != -1)
                    {
                        c = v.label;
                        v.TexCoordinate.x = bu[c];
                        v.TexCoordinate.y = bv[c];
                        v.TexCoordinate.z = 0.0;
                    }
                }
            }
            Console.Write("Stretchmin done\n");
            Copy(opt);
            Console.Write("opt assigned\n");
        }


        #region Subdivide

        //Loop subdivision only works on triangle meshes
        //Loop subdivision only works on 2-manifold meshes
        public Mesh SubdivideLoop()
        {
            int i, j, k;
            Vector3 nv;
            Mesh nm = new Mesh();

            for (i = 0; i < Vertices.Count; i++)
            {
                nv = new Vector3(0.0, 0.0, 0.0);
                if (Vertices[i].IsBoundary)
                {
                    //find two boundary edges
                    // each opposite end contributes 1/8
                    // while this vertex contributes 6/8
                    for (j = 0; j < Vertices[i].Edges.Count; j++)
                    {
                        Edge e = Vertices[i].Edges[j];
                        if (e.IsBoundary)
                        {
                            nv += 1.0 / 8.0 * e.Vertices[0].Coordinate;
                            nv += 1.0 / 8.0 * e.Vertices[1].Coordinate;
                        }
                    }
                    //this vertex has been added with weight 2*1.0/8.0
                    //still need 4.0/8.0
                    nv += 0.5 * Vertices[i].Coordinate;
                }
                else
                {
                    UnicityList<Vertex> ring = ComputeExclusiveRing(Vertices[i], 1);
                    //the valence
                    int n = Vertices[i].Edges.Count; //==ring.getCount()
                    //Loop's method
                    //double num=(3.0+2.0*System.Math.Cos(2.0*System.Math.PI/n));
                    //double beta=1.0/n*(5.0/8.0-num*num/64.0);
                    //Warren and Weimer's method			                   
                    double beta = 3.0 / (n * (n + 2.0));
                    nv = Vertices[i].Coordinate * (1 - n * beta);
                    for (j = 0; j < n; j++)
                    {
                        nv += beta * ring[j].Coordinate;
                    }
                }
                Vertices[i].label = nm.AddVertex(nv);
            }
            for (i = 0; i < Edges.Count; i++)
            {
                nv = new Vector3();
                Edge e = Edges[i];
                if (e.IsBoundary)
                {
                    nv = 0.5 * e.Vertices[0].Coordinate;
                    nv += 0.5 * e.Vertices[1].Coordinate;
                }
                else
                {
                    nv = 3.0 / 8.0 * e.Vertices[0].Coordinate;
                    nv += 3.0 / 8.0 * e.Vertices[1].Coordinate;
                    //for each face find the vertex
                    //not on the edge
                    for (k = 0; k < 2; k++)
                    {
                        nv += 1.0 / 8.0 * e.Faces[k].OffEdge(e).Coordinate;
                    }
                }
                e.label = nm.AddVertex(nv);
            }

            //vertices are created, now add the triangles
            for (i = 0; i < Faces.Count; i++)
            {
                Face f = Faces[i];
                //center triangle
                nm.AddTriangle(f.Edges[0].label, f.Edges[1].label, f.Edges[2].label);
                //triangles connected to existing vertices
                for (j = 0; j < 3; j++)
                {
                    for (k = j + 1; k < 3; k++)
                    {
                        Vertex v = null;
                        //find common vertex
                        if (f.Edges[j].Vertices[0] == f.Edges[k].Vertices[0])
                            v = f.Edges[j].Vertices[0];
                        if (f.Edges[j].Vertices[0] == f.Edges[k].Vertices[1])
                            v = f.Edges[j].Vertices[0];
                        if (f.Edges[j].Vertices[1] == f.Edges[k].Vertices[1])
                            v = f.Edges[j].Vertices[1];
                        if (f.Edges[j].Vertices[1] == f.Edges[k].Vertices[0])
                            v = f.Edges[j].Vertices[1];
                        nm.AddTriangle(f.Edges[j].label, v.label, f.Edges[k].label);
                    }
                }
            }

            return nm;
        }

        //this is an interpolating subdivision scheme
        // so existing vertices are left unaltered
        public Mesh SubdivideButterfly()
        {
            int i, j, k;
            Vector3 nv;
            Mesh nm = new Mesh();

            //add existing vertices with no changes
            for (i = 0; i < Vertices.Count; i++)
            {
                Vertices[i].label = nm.AddVertex(Vertices[i].Coordinate);
            }
            //create new vertices for each edge
            for (i = 0; i < Edges.Count; i++)
            {
                nv = new Vector3();
                Edge e = Edges[i];
                if (e.IsBoundary)
                {
                    nv = 9.0 / 16.0 * e.Vertices[0].Coordinate;
                    nv += 9.0 / 16.0 * e.Vertices[1].Coordinate;
                    //find adjacent boundary edges
                    for (k = 0; k < 2; k++)
                    {
                        for (j = 0; j < e.Vertices[k].Edges.Count; j++)
                        {
                            Edge eb = e.Vertices[k].Edges[j];
                            if ((eb != e) && (eb.IsBoundary))
                            {
                                if (eb.Vertices[0] == e.Vertices[k])
                                    nv -= 1.0 / 16.0 * eb.Vertices[1].Coordinate;
                                else
                                    nv -= 1.0 / 16.0 * eb.Vertices[0].Coordinate;
                            }
                        }
                    }
                }
                else
                {
                    //number of irregular settings encountered
                    // to compute average
                    int irreg = 0;
                    //check for regular setting
                    if ((e.Vertices[0].Edges.Count == 6)
                        && (e.Vertices[1].Edges.Count == 6))
                    {
                        nv += 0.5 * e.Vertices[0].Coordinate;
                        nv += 0.5 * e.Vertices[1].Coordinate;
                        //for each face find the vertex
                        //not on the edge
                        for (k = 0; k < 2; k++)
                        {
                            nv += 1.0 / 8.0 * e.Faces[k].OffEdge(e).Coordinate;
                        }
                        //Now find the other faces adjacent to the
                        // edge faces, and find the point not on the
                        // edge faces (4 of them).
                        for (k = 0; k < 2; k++)
                        {
                            for (j = 0; j < 3; j++)
                            {
                                if (e.Faces[k].Edges[j] != e)
                                {
                                    Face adj = null;
                                    if (e.Faces[k] == e.Faces[k].Edges[j].Faces[0])
                                        adj = e.Faces[k].Edges[j].Faces[1];
                                    else
                                        adj = e.Faces[k].Edges[j].Faces[0];
                                    if (adj != null)
                                        nv -= 1.0 / 16.0 * adj.OffEdge(e.Faces[k].Edges[j]).Coordinate;
                                    else
                                    {
                                        //use a reflection
                                    }
                                }
                            }
                        }

                        //since we will divide set irreg to 
                        // something not zero
                        irreg = 1;
                    }
                    //and now the irregular setting
                    for (k = 0; k < 2; k++)
                    {
                        if (e.Vertices[k].Edges.Count != 6)
                        {
                            double sumwi = 0.0;
                            UnicityList<Vertex> ring = ComputeExclusiveRing(e.Vertices[k], 1);
                            ring = OrderRing(e.Vertices[k], e.Vertices[1 - k], ring);
                            if (e.Vertices[k].Edges.Count == 3)
                            {
                                nv += 5.0 / 12.0 * ring[0].Coordinate; sumwi += 5.0 / 12.0;
                                nv -= 1.0 / 12.0 * ring[1].Coordinate; sumwi += -1.0 / 12.0;
                                nv -= 1.0 / 12.0 * ring[2].Coordinate; sumwi += -1.0 / 12.0;
                            }
                            else
                                if (e.Vertices[k].Edges.Count == 4)
                                {
                                    nv += 3.0 / 8.0 * ring[0].Coordinate; sumwi += 3.0 / 8.0;
                                    //ring[1] weight is zero
                                    nv -= 1.0 / 8.0 * ring[2].Coordinate; sumwi += -1.0 / 8.0;
                                    //ring[3] weight is zero
                                }
                                else
                                {
                                    int n = ring.Count;
                                    for (j = 0; j < ring.Count; j++)
                                    {
                                        double wi = 0.25 + System.Math.Cos(2.0 * System.Math.PI * j / n) +
                                            0.5 * System.Math.Cos(4.0 * System.Math.PI * j / n);
                                        wi /= n;
                                        sumwi += wi;
                                        nv += wi * ring[j].Coordinate;
                                    }
                                }
                            nv += (1.0 - sumwi) * e.Vertices[k].Coordinate;
                            irreg++;
                        }
                    }
                    nv /= irreg;
                }
                e.label = nm.AddVertex(nv);
            }

            //vertices are created, now add the triangles
            for (i = 0; i < Faces.Count; i++)
            {
                Face f = Faces[i];
                //center triangle
                nm.AddTriangle(f.Edges[0].label, f.Edges[1].label, f.Edges[2].label);
                //triangles connected to existing vertices
                for (j = 0; j < 3; j++)
                {
                    for (k = j + 1; k < 3; k++)
                    {
                        Vertex v = null;
                        //find common vertex
                        if (f.Edges[j].Vertices[0] == f.Edges[k].Vertices[0])
                            v = f.Edges[j].Vertices[0];
                        if (f.Edges[j].Vertices[0] == f.Edges[k].Vertices[1])
                            v = f.Edges[j].Vertices[0];
                        if (f.Edges[j].Vertices[1] == f.Edges[k].Vertices[1])
                            v = f.Edges[j].Vertices[1];
                        if (f.Edges[j].Vertices[1] == f.Edges[k].Vertices[0])
                            v = f.Edges[j].Vertices[1];
                        nm.AddTriangle(f.Edges[j].label, v.label, f.Edges[k].label);
                    }
                }
            }

            return nm;
        }

        #endregion

        //boundary specifies whether the boundary edge must be subdivided
        //this must only happen every second iteration
        public Mesh SubdivideSqrt3(bool boundary)
        {
            int i, j, k;
            int v1, v2, v3, v4, v5;
            Mesh nm = new Mesh();
            //edge flipping happens as a second pass
            //so we have two meshes
            Mesh nm2 = new Mesh();
            Vector3 nv;
            //smooth exisitng vertices
            for (i = 0; i < Vertices.Count; i++)
            {
                UnicityList<Vertex> ring;
                nv = new Vector3();
                if (Vertices[i].IsBoundary)
                {
                    if (boundary)
                    {
                        //find two boundary edges
                        // each opposite end contributes 4/27
                        // while this vertex contributes 19/27
                        for (j = 0; j < Vertices[i].Edges.Count; j++)
                        {
                            Edge e = Vertices[i].Edges[j];
                            if (e.IsBoundary)
                            {
                                nv += 4.0 / 27.0 * e.Vertices[0].Coordinate;
                                nv += 4.0 / 27.0 * e.Vertices[1].Coordinate;
                            }
                        }
                        //this vertex has been added with weight 2*4.0/27.0
                        //still need 11.0/27.0
                        nv += 11.0 / 27.0 * Vertices[i].Coordinate;
                    }
                    else
                    {
                        nv = (Vector3)Vertices[i].Coordinate.Clone();
                    }
                }
                else
                {
                    int n = Vertices[i].Edges.Count;
                    double alpha = 4.0 - 2.0 * System.Math.Cos(2.0 * System.Math.PI / n);
                    alpha /= 9.0;
                    nv = (1.0 - alpha) * Vertices[i].Coordinate;
                    ring = ComputeExclusiveRing(Vertices[i], 1);
                    for (j = 0; j < ring.Count; j++)
                    {
                        nv += alpha / n * ring[j].Coordinate;
                    }
                }
                Vertices[i].label = nm.AddVertex(nv);
                //label this as : not a new vertex
                nm.Vertices[Vertices[i].label].label = 0;
            }
            //add a new vertex for each face
            for (i = 0; i < Faces.Count; i++)
            {
                nv = new Vector3();
                if ((Faces[i].IsBoundary) && (boundary))
                {
                    Edge e = null;
                    //there should be precisely one boundary edge
                    //in the first pass any triangle with two boundary
                    // edges is subdivided so that each new triangle
                    // has only one boundary edge. Since this code
                    // should only be invoked every second iteration
                    // only one bounadry edge should be present.
                    for (j = 0; j < Faces[i].Edges.Count; j++)
                    {
                        if (Faces[i].Edges[j].IsBoundary)
                            e = Faces[i].Edges[j];
                    }
                    //split the boundary edge
                    //find two boundary edges
                    nv = new Vector3();
                    for (k = 0; k < e.Vertices[0].Edges.Count; k++)
                    {
                        Edge e2 = e.Vertices[0].Edges[k];
                        if (e2.IsBoundary && (e2 != e))
                        {
                            if (e2.Vertices[0] != e.Vertices[0])
                            {
                                nv = 1.0 / 27.0 * (e2.Vertices[0].Coordinate +
                                             16.0 * e.Vertices[0].Coordinate +
                                             10.0 * e.Vertices[1].Coordinate);
                            }
                            else
                            {
                                nv = 1.0 / 27.0 * (e2.Vertices[1].Coordinate +
                                             16.0 * e.Vertices[0].Coordinate +
                                             10.0 * e.Vertices[1].Coordinate);
                            }
                        }
                    }
                    //nv=e.ends[0].v*2.0/3.0+e.ends[1].v*1.0/3.0;
                    v4 = nm.AddVertex(nv);

                    nv = new Vector3();
                    for (k = 0; k < e.Vertices[1].Edges.Count; k++)
                    {
                        Edge e2 = e.Vertices[1].Edges[k];
                        if (e2.IsBoundary && (e2 != e))
                        {
                            if (e2.Vertices[0] != e.Vertices[1])
                            {
                                nv = 1.0 / 27.0 * (e2.Vertices[0].Coordinate +
                                             16.0 * e.Vertices[1].Coordinate +
                                             10.0 * e.Vertices[0].Coordinate);
                            }
                            else
                            {
                                nv = 1.0 / 27.0 * (e2.Vertices[1].Coordinate +
                                             16.0 * e.Vertices[1].Coordinate +
                                             10.0 * e.Vertices[0].Coordinate);
                            }
                        }
                    }
                    //nv=e.ends[0].v*1.0/3.0+e.ends[1].v*2.0/3.0;
                    v5 = nm.AddVertex(nv);
                    //since vertices are added sequntially we have v5=v4+1
                    Faces[i].label = v4;
                    nm.Vertices[v4].label = 1;
                    nm.Vertices[v5].label = 1;
                }
                else
                {
                    nv = Faces[i].Vertices[0].Coordinate / 3.0;
                    nv += Faces[i].Vertices[1].Coordinate / 3.0;
                    nv += Faces[i].Vertices[2].Coordinate / 3.0;
                    Faces[i].label = nm.AddVertex(nv);
                    //label this as : a new vertex
                    nm.Vertices[Faces[i].label].label = 1;
                }
            }
            //now create the new triangles
            //edge flips are done later
            for (i = 0; i < Faces.Count; i++)
            {
                if ((Faces[i].IsBoundary) && (boundary))
                {
                    //boundary edges are only subdivided
                    // every second iteration and are not subject
                    // to edge flips.
                    Edge e = null;
                    //there should be precisely one boundary edge
                    //in the first pass any triangle with two boundary
                    // edges is subdivided so that each new triangle
                    // has only one boundary edge. Since this code
                    // should only be invoked every second iteration
                    // only one bounadry edge should be present.
                    for (j = 0; j < Faces[i].Edges.Count; j++)
                    {
                        if (Faces[i].Edges[j].IsBoundary)
                            e = Faces[i].Edges[j];
                    }
                    v1 = Faces[i].OffEdge(e).label;
                    v2 = e.Vertices[0].label;
                    v3 = e.Vertices[1].label;
                    v4 = Faces[i].label;
                    v5 = v4 + 1;
                    nm.AddTriangle(v1, v2, v4);
                    nm.AddTriangle(v1, v4, v5);
                    nm.AddTriangle(v1, v5, v3);
                }
                else
                {
                    v4 = Faces[i].label;
                    v1 = Faces[i].Vertices[0].label;
                    v2 = Faces[i].Vertices[1].label;
                    v3 = Faces[i].Vertices[2].label;
                    nm.AddTriangle(v1, v2, v4);
                    nm.AddTriangle(v2, v3, v4);
                    nm.AddTriangle(v3, v1, v4);
                }
            }

            //now do edge flipping, create a new mesh with the
            // same vertices
            for (i = 0; i < nm.Vertices.Count; i++)
            {
                //add the vertex and store the index
                nm.Vertices[i].IndexInMesh = nm2.AddVertex(nm.Vertices[i].Coordinate);
            }
            //now add the triangles in such a way that the edges
            // are flipped
            for (i = 0; i < nm.Faces.Count; i++)
            {
                //not added yet
                nm.Faces[i].label = 0;
            }
            for (i = 0; i < nm.Edges.Count; i++)
            {
                //only non boundary edges need to be flipped
                Edge e = nm.Edges[i];
                if (!e.IsBoundary)
                {
                    //only flip if both ends of the edge are "old vertices"
                    if ((e.Vertices[0].label == 0) &&
                        (e.Vertices[1].label == 0))
                    {
                        //and if oppposite ends are "new vertices"
                        if ((e.Faces[0].OffEdge(e).label == 1)
                            && (e.Faces[1].OffEdge(e).label == 1))
                        {
                            v1 = e.Vertices[0].IndexInMesh;
                            v2 = e.Vertices[1].IndexInMesh;
                            v3 = e.Faces[0].OffEdge(e).IndexInMesh;
                            v4 = e.Faces[1].OffEdge(e).IndexInMesh;
                            nm2.AddTriangle(v1, v3, v4);
                            nm2.AddTriangle(v3, v2, v4);
                            //these triangles are added
                            e.Faces[0].label = 1;
                            e.Faces[1].label = 1;
                        }
                    }
                }
            }
            for (i = 0; i < nm.Faces.Count; i++)
            {
                if (nm.Faces[i].label == 0)
                {
                    v1 = nm.Faces[i].Vertices[0].IndexInMesh;
                    v2 = nm.Faces[i].Vertices[1].IndexInMesh;
                    v3 = nm.Faces[i].Vertices[2].IndexInMesh;
                    nm2.AddTriangle(v1, v2, v3);
                }
            }
            return nm2;
        }

        #region AlignAdjacentFaces

        public void AlignAdjacentFaces(UnicityList<Face> queue)
        {
            int i;
            Face f;

            while (queue.Count > 0)
            {
                f = queue[0];
                queue.Remove(0);
                if (f.label == 0) return;
                for (i = 0; i < f.Edges.Count; i++)
                {
                    Edge e = f.Edges[i];
                    Face adj;
                    if (!e.IsBoundary)
                    {
                        if (e.Faces[0] == f) adj = e.Faces[1];
                        else adj = e.Faces[0];
                        //might cause acute angles to face inward!
                        if (adj.Normal.Dot(f.Normal) < 0.0 && adj.label == 0)
                        {
                            //flip adj to have same orientation
                            // as this face;
                            Vertex tmp;
                            tmp = adj.Vertices[0];
                            adj.Vertices[0] = adj.Vertices[2];
                            adj.Vertices[2] = tmp;
                            adj.Normal = -adj.Normal;
                        }
                        if (adj.label == 0)
                        {
                            adj.label = 1;
                            queue.Add(adj);
                        }
                    }
                }
            }
        }

        #endregion

        #region AlignNormals

        //try to fix problems in alignment
        // of normals caused by incorrect
        // triangle vertex orders
        //Only works on triangles
        public void AlignNormals()
        {
            int i;
            ComputeFaceNormals();
            for (i = 0; i < Faces.Count; i++)
                Faces[i].label = 0;

            Faces[0].label = 1;
            UnicityList<Face> queue = new UnicityList<Face>();
            queue.Add(Faces[0]);
            AlignAdjacentFaces(queue);
            ComputeVertexNormals();
        }

        #endregion

        public void SaveTexInNormal()
        {
            for (int i = 0; i < Vertices.Count; i++)
                Vertices[i].Normal = Vertices[i].TexCoordinate;
        }

        #region Wavefront

        //public void WriteOFF(string filename)
        //{
        //    int nv, nf;
        //    int i, v1, v2, v3;
        //    StreamWriter f = new StreamWriter(filename);
        //    if (f == null) return;

        //    f.Write("OFF\n");
        //    nv = vertices.getCount();
        //    nf = faces.getCount();
        //    f.Write("{0} {1} 0\n", nv, nf);
        //    for (i = 0; i < nv; i++)
        //    {
        //        f.Write("{0} {1} {2}\n", vertices[i].v.x, vertices[i].v.y, vertices[i].v.z);
        //        vertices[i].label = i;
        //    }
        //    for (i = 0; i < nf; i++)
        //    {
        //        v1 = faces[i].vertices[0].label;
        //        v2 = faces[i].vertices[1].label;
        //        v3 = faces[i].vertices[2].label;
        //        f.Write("3 {0} {1} {2}\n", v1, v2, v3);
        //    }
        //    f.Close();
        //}

        ////wavefront obj
        //public void WriteOBJ(string filename)
        //{
        //    int nv, nf;
        //    int i, v1, v2, v3;
        //    StreamWriter f = new StreamWriter(filename);
        //    if (f == null) return;

        //    nv = vertices.getCount();
        //    nf = faces.getCount();
        //    for (i = 0; i < nv; i++)
        //    {
        //        f.Write("v {0} {1} {2}\n", vertices[i].v.x, vertices[i].v.y, vertices[i].v.z);
        //        vertices[i].label = i + 1;
        //    }
        //    for (i = 0; i < nf; i++)
        //    {
        //        v1 = faces[i].vertices[0].label;
        //        v2 = faces[i].vertices[1].label;
        //        v3 = faces[i].vertices[2].label;
        //        f.Write("f {0} {1} {2}\n", v1, v2, v3);
        //    }
        //    f.Close();
        //}

        //public void WriteParamOBJ(string filename)
        //{
        //    int nv, nf;
        //    int i, v1, v2, v3;
        //    StreamWriter f = new StreamWriter(filename);
        //    if (f == null) return;

        //    nv = vertices.getCount();
        //    nf = faces.getCount();
        //    for (i = 0; i < nv; i++)
        //    {
        //        f.Write("v {0} {1} {2}\n", vertices[i].tex.x, vertices[i].tex.y, vertices[i].tex.z);
        //        vertices[i].label = i + 1;
        //    }
        //    for (i = 0; i < nf; i++)
        //    {
        //        v1 = faces[i].vertices[0].label;
        //        v2 = faces[i].vertices[1].label;
        //        v3 = faces[i].vertices[2].label;
        //        f.Write("f {0} {1} {2}\n", v1, v2, v3);
        //    }
        //    f.Close();
        //}

        //public void WriteParamOFF(string filename)
        //{
        //    int nv, nf;
        //    int i, v1, v2, v3;
        //    StreamWriter f = new StreamWriter(filename);
        //    if (f == null) return;

        //    f.Write("OFF\n");
        //    nv = vertices.getCount();
        //    nf = faces.getCount();
        //    f.Write("{0} {1} 0\n", nv, nf);
        //    for (i = 0; i < nv; i++)
        //    {
        //        f.Write("{0} {1} {2}\n", vertices[i].tex.x, vertices[i].tex.y, vertices[i].tex.z);
        //        vertices[i].label = i;
        //    }
        //    for (i = 0; i < nf; i++)
        //    {
        //        v1 = faces[i].vertices[0].label;
        //        v2 = faces[i].vertices[1].label;
        //        v3 = faces[i].vertices[2].label;
        //        f.Write("3 {0} {1} {2}\n", v1, v2, v3);
        //    }
        //    f.Close();
        //}

        #endregion

        // Size is the recommended size, but it may have to be larger to avoid
        // vertices being placed in the same cell
        public Vector3[,] CreateGeometryImage(int size, bool resize)
        {
            return CreateGeometryImage(size, resize, GEOMETRY_IMAGE);
        }

        public Vector3[,] CreateGeometryImage(int size, bool resize, int type)
        {
            int[,] occupied;
            Vector3[,] gim = null;
            int i, j, x, y, newsize, closest;
            bool restart = true;
            double minx, maxx, miny, maxy, tx, ty;
            double alpha, beta, gamma;
            double basearea, min;
            Face f;

            newsize = size;
            while (restart)
            {
                restart = false;
                size = newsize;
                //no vertices have been placed
                occupied = new int[size, size];
                for (y = 0; y < size; y++)
                {
                    for (x = 0; x < size; x++)
                    {
                        occupied[y, x] = -1;
                    }
                }

                //place the vertices in cells				
                for (i = 0; i < Vertices.Count; i++)
                {
                    x = (int)(Vertices[i].TexCoordinate.x * (size - 1) + 0.4);
                    y = (int)(Vertices[i].TexCoordinate.y * (size - 1) + 0.4);
                    //if not occupied, then use the vertex
                    if (occupied[y, x] == -1)
                    {
                        occupied[y, x] = i;
                    }
                    else
                    {
                        //occupied cell, either don't use the
                        //   current vertex, or resize and start again
                        if ((resize) && (size < 1024))
                        {
                            newsize = size * 2;
                            restart = true;
                            break;
                        }
                    }
                }

                //if there is no need to restart, then search through
                // cells for unoccupied cells
                if (!restart)
                {
                    gim = new Vector3[size, size];
                    for (y = 0; y < size; y++)
                    {
                        UnicityList<Face> linecandidates = new UnicityList<Face>();
                        ty = ((y) / (double)(size - 1));
                        if (y == 0) ty = 0.0;
                        if (y == size - 1) ty = 1.0;

                        for (i = 0; i < Faces.Count; i++)
                        {
                            //compute the bounding box of the triangle
                            //minx=maxx=faces[i].vertices[0].tex.x;
                            miny = maxy = Faces[i].Vertices[0].TexCoordinate.y;
                            for (j = 1; j < 3; j++)
                            {
                                //if (faces[i].vertices[j].tex.x>maxx)
                                //	maxx=faces[i].vertices[j].tex.x;
                                //if (faces[i].vertices[j].tex.x<minx)
                                //	minx=faces[i].vertices[j].tex.x;
                                if (Faces[i].Vertices[j].TexCoordinate.y > maxy)
                                    maxy = Faces[i].Vertices[j].TexCoordinate.y;
                                if (Faces[i].Vertices[j].TexCoordinate.y < miny)
                                    miny = Faces[i].Vertices[j].TexCoordinate.y;
                            }
                            //check to see if it is in the bounding box
                            if (/*(minx<=tx)&&(tx<=maxx)&&*/(miny <= ty) && (ty <= maxy))
                            {
                                linecandidates.Add(Faces[i]);
                            }
                        }

                        for (x = 0; x < size; x++)
                        {
                            if (occupied[y, x] != -1)
                            {
                                if (type == NORMAL_MAP)
                                {
                                    gim[y, x] = Vertices[occupied[y, x]].Normal;
                                }
                                else
                                {
                                    gim[y, x] = Vertices[occupied[y, x]].Coordinate;
                                }
                            }
                            else
                                if (occupied[y, x] == -1)
                                {
                                    //take the unnocupied cell	
                                    // and compute in which triangle the
                                    // centre of the cell lies
                                    UnicityList<Face> candidates = new UnicityList<Face>();
                                    //tx=((x+0.5)/(double)size);
                                    //ty=((y+0.5)/(double)size);
                                    tx = ((x) / (double)(size - 1));
                                    //ty=((y)/(double)(size-1));
                                    // handle borders specially
                                    if (x == 0) tx = 0.0;
                                    if (x == size - 1) tx = 1.0;
                                    //if (y==0) ty=0.0;
                                    //if (y==size-1) ty=1.0;
                                    Vector3 v = new Vector3();
                                    v.x = tx; v.y = ty; v.z = 0.0;
                                    for (i = 0; i < linecandidates.Count; i++)
                                    {
                                        //compute the bounding box of the triangle
                                        minx = maxx = linecandidates[i].Vertices[0].TexCoordinate.x;
                                        //miny=maxy=faces[i].vertices[0].tex.y;
                                        for (j = 1; j < 3; j++)
                                        {
                                            if (linecandidates[i].Vertices[j].TexCoordinate.x > maxx)
                                                maxx = linecandidates[i].Vertices[j].TexCoordinate.x;
                                            if (linecandidates[i].Vertices[j].TexCoordinate.x < minx)
                                                minx = linecandidates[i].Vertices[j].TexCoordinate.x;
                                            //if (faces[i].vertices[j].tex.y>maxy)
                                            //	maxy=faces[i].vertices[j].tex.y;
                                            //if (faces[i].vertices[j].tex.y<miny)
                                            //	miny=faces[i].vertices[j].tex.y;
                                        }
                                        //check to see if it is in the bounding box
                                        if ((minx <= tx) && (tx <= maxx)/*&&(miny<=ty)&&(ty<=maxy)*/)
                                        {
                                            //candidates.add(faces[i]);
                                            candidates.Add(linecandidates[i]);
                                        }
                                    }
                                    closest = -1;
                                    min = 1e10;
                                    Vector3 v1, v2, v3;
                                    for (i = 0; i < candidates.Count; i++)
                                    {
                                        //compute barycentric coordinates
                                        // and use them to determine the 
                                        // euclidean 3d-space position
                                        //the point is inside the triangle if
                                        // the sum of the barycentric coordinates is 1
                                        // and all coordinates are greater than 0.
                                        f = candidates[i];
                                        v1 = f.Vertices[0].TexCoordinate;
                                        v2 = f.Vertices[1].TexCoordinate;
                                        v3 = f.Vertices[2].TexCoordinate;
                                        basearea = area(v1, v2, v3);
                                        alpha = area(v, v2, v3) / basearea;
                                        beta = area(v1, v, v3) / basearea;
                                        gamma = area(v1, v2, v) / basearea;
                                        //alpha, beta and gamma are all positive by definition	
                                        //test if they add to 1.0
                                        //if (System.Math.Abs(alpha+beta+gamma-1.0)<1e-7) {
                                        //	v=alpha*f.vertices[0].v+
                                        //		beta*f.vertices[1].v+	
                                        //		gamma*f.vertices[2].v;	
                                        //} else {
                                        if (System.Math.Abs(alpha + beta + gamma - 1.0) < min)
                                        {
                                            closest = i;
                                            min = System.Math.Abs(alpha + beta + gamma - 1.0);
                                        }
                                        //}
                                    }
                                    f = candidates[closest];
                                    v1 = f.Vertices[0].TexCoordinate;
                                    v2 = f.Vertices[1].TexCoordinate;
                                    v3 = f.Vertices[2].TexCoordinate;
                                    basearea = area(v1, v2, v3);
                                    alpha = area(v, v2, v3) / basearea;
                                    beta = area(v1, v, v3) / basearea;
                                    gamma = area(v1, v2, v) / basearea;
                                    if (type == NORMAL_MAP)
                                    {
                                        v = alpha * f.Vertices[0].Normal +
                                            beta * f.Vertices[1].Normal +
                                            gamma * f.Vertices[2].Normal;
                                    }
                                    else
                                    {
                                        v = alpha * f.Vertices[0].Coordinate +
                                            beta * f.Vertices[1].Coordinate +
                                            gamma * f.Vertices[2].Coordinate;
                                    }
                                    gim[y, x] = v;
                                }
                        }
                    }
                }
            }

            return gim;
        }

        #region RIB

        //public void WriteGeometryImageRIB(string filename, string tiffname, int size)
        //{
        //    int i;
        //    StreamWriter f = new StreamWriter(filename);
        //    Vector3 v1, v2, v3, v4;

        //    if (f == null) return;
        //    f.Write("Format {0} {1} 1\n", size, size);
        //    f.Write("Clipping 0.1 100.0\n");
        //    f.Write("PixelSamples 2 2\n");
        //    f.Write("Sides 2\n");
        //    f.Write("ShadingRate 0.5\n");
        //    f.Write("Quantize \"rgba\" 0 0 0 0\n");
        //    f.Write("Display \"{0}\" \"file\" \"rgba\"\n", tiffname);
        //    f.Write("Scale {0} {1} 1\n", 2.0 - 1.0 / size, 2.0 - 1.0 / size);
        //    f.Write("Translate 0 0 1\n");
        //    f.Write("WorldBegin\n");
        //    f.Write("\tSurface \"constant\"\n");

        //    for (i = 0; i < vertices.getCount(); i++)
        //    {
        //        vertices[i].tex.x -= 0.5;
        //        vertices[i].tex.y -= 0.5;
        //        vertices[i].label = i;
        //    }
        //    f.Write("\t PointsPolygons [\n");
        //    for (i = 0; i < faces.getCount(); i++)
        //    {
        //        f.Write("\t\t 3\n");
        //    }
        //    f.Write("\t ] [\n");
        //    for (i = 0; i < faces.getCount(); i++)
        //    {
        //        f.Write("\t\t {0} {1} {2}\n", faces[i].vertices[0].label,
        //                faces[i].vertices[1].label,
        //                faces[i].vertices[2].label);
        //    }
        //    f.Write("\t ]\n");
        //    f.Write("\t \"P\" [\n");
        //    for (i = 0; i < vertices.getCount(); i++)
        //    {
        //        f.Write("\t\t {0} {1} {2} \n",
        //            vertices[i].tex.x,
        //            vertices[i].tex.y,
        //            vertices[i].tex.z);
        //    }
        //    f.Write("\t ]\n");
        //    f.Write("\t \"Cs\" [\n");
        //    for (i = 0; i < vertices.getCount(); i++)
        //    {
        //        f.Write("\t\t {0} {1} {2} \n",
        //            vertices[i].v.x,
        //            vertices[i].v.y,
        //            vertices[i].v.z);
        //    }
        //    f.Write("\t ]\n");
        //    f.Write("\tTranslate 0.0 0.0 0.5\n");
        //    for (i = 0; i < edges.getCount(); i++)
        //    {
        //        if (edges[i].boundary())
        //        {
        //            v1 = edges[i].ends[0].tex * 0.5;
        //            v2 = edges[i].ends[0].tex * 1.5;
        //            v3 = edges[i].ends[1].tex * 1.5;
        //            v4 = edges[i].ends[1].tex * 0.5;
        //            f.Write("\t Polygon \"Cs\" [ ");
        //            f.Write("{0} {1} {2} ",
        //                edges[i].ends[0].v.x,
        //                edges[i].ends[0].v.y,
        //                edges[i].ends[0].v.z);
        //            f.Write("{0} {1} {2} ",
        //                edges[i].ends[0].v.x,
        //                edges[i].ends[0].v.y,
        //                edges[i].ends[0].v.z);
        //            f.Write("{0} {1} {2} ",
        //                edges[i].ends[1].v.x,
        //                edges[i].ends[1].v.y,
        //                edges[i].ends[1].v.z);
        //            f.Write("{0} {1} {2} ",
        //                edges[i].ends[1].v.x,
        //                edges[i].ends[1].v.y,
        //                edges[i].ends[1].v.z);
        //            f.Write("]\n\t\t \"P\" [ ");
        //            f.Write("{0} {1} {2} ", v1.x, v1.y, v1.z);
        //            f.Write("{0} {1} {2} ", v2.x, v2.y, v2.z);
        //            f.Write("{0} {1} {2} ", v3.x, v3.y, v3.z);
        //            f.Write("{0} {1} {2} ", v4.x, v4.y, v4.z);
        //            f.Write("]\n");
        //        }
        //    }
        //    f.Write("WorldEnd\n");
        //    f.Close();

        //    for (i = 0; i < vertices.getCount(); i++)
        //    {
        //        vertices[i].tex.x += 0.5;
        //        vertices[i].tex.y += 0.5;
        //    }
        //}

        //public void WriteParamMetaPost(string filename)
        //{
        //    int i;
        //    StreamWriter f = new StreamWriter(filename);
        //    if (f == null) return;
        //    Vector3 v1, v2;

        //    f.Write("u:=5cm;\n");
        //    f.Write("beginfig(1);\n");
        //    for (i = 0; i < edges.getCount(); i++)
        //    {
        //        v1 = edges[i].ends[0].tex;
        //        v2 = edges[i].ends[1].tex;
        //        f.Write("\tdraw ({0},{1})--({0},{1});\n", v1.x, v1.y, v2.x, v2.y);
        //    }
        //    f.Write("endfig;\n");
        //    f.Write("end;\n");
        //    f.Close();
        //}

        //public void WriteNormalMapRIB(string filename, string tiffname, int size)
        //{
        //    int i;

        //    StreamWriter f = new StreamWriter(filename);
        //    Vector3 v1, v2, v3, v4;
        //    if (f == null) return;
        //    f.Write("Format {0} {1} 1\n", size, size);
        //    f.Write("Clipping 0.1 100.0\n");
        //    f.Write("PixelSamples 2 2\n");
        //    f.Write("Sides 2\n");
        //    f.Write("ShadingRate 0.5\n");
        //    f.Write("Quantize \"rgba\" 0 0 0 0\n");
        //    f.Write("Display \"{0}\" \"file\" \"rgba\"\n", tiffname);
        //    f.Write("Scale {0} {1} 1\n", 2.0 - 1.0 / size, 2.0 - 1.0 / size);
        //    f.Write("Translate 0 0 1\n");
        //    f.Write("WorldBegin\n");
        //    f.Write("\tSurface \"constant\"\n");
        //    for (i = 0; i < vertices.getCount(); i++)
        //    {
        //        vertices[i].tex.x -= 0.5;
        //        vertices[i].tex.y -= 0.5;
        //        vertices[i].label = i;
        //    }

        //    f.Write("\t PointsPolygons [\n");
        //    for (i = 0; i < faces.getCount(); i++)
        //    {
        //        f.Write("\t\t 3\n");
        //    }
        //    f.Write("\t ] [\n");
        //    for (i = 0; i < faces.getCount(); i++)
        //    {
        //        f.Write("\t\t {0} {1} {2}\n", faces[i].vertices[0].label,
        //                faces[i].vertices[1].label,
        //                faces[i].vertices[2].label);
        //    }
        //    f.Write("\t ]\n");
        //    f.Write("\t \"P\" [\n");
        //    for (i = 0; i < vertices.getCount(); i++)
        //    {
        //        f.Write("\t\t {0} {1} {2} \n",
        //            vertices[i].tex.x,
        //            vertices[i].tex.y,
        //            vertices[i].tex.z);
        //    }
        //    f.Write("\t ]\n");
        //    f.Write("\t \"Cs\" [\n");
        //    for (i = 0; i < vertices.getCount(); i++)
        //    {
        //        f.Write("\t\t {0} {1} {2} \n",
        //            vertices[i].normal.x,
        //            vertices[i].normal.y,
        //            vertices[i].normal.z);
        //    }
        //    f.Write("\t ]\n");
        //    f.Write("\tTranslate 0.0 0.0 0.5\n");
        //    for (i = 0; i < edges.getCount(); i++)
        //    {
        //        if (edges[i].boundary())
        //        {
        //            v1 = edges[i].ends[0].tex * 0.5;
        //            v2 = edges[i].ends[0].tex * 1.5;
        //            v3 = edges[i].ends[1].tex * 1.5;
        //            v4 = edges[i].ends[1].tex * 0.5;
        //            f.Write("\t Polygon \"Cs\" [ ");
        //            f.Write("{0} {1} {2} ",
        //                edges[i].ends[0].normal.x,
        //                edges[i].ends[0].normal.y,
        //                edges[i].ends[0].normal.z);
        //            f.Write("{0} {1} {2} ",
        //                edges[i].ends[0].normal.x,
        //                edges[i].ends[0].normal.y,
        //                edges[i].ends[0].normal.z);
        //            f.Write("{0} {1} {2} ",
        //                edges[i].ends[1].normal.x,
        //                edges[i].ends[1].normal.y,
        //                edges[i].ends[1].normal.z);
        //            f.Write("{0} {1} {2} ",
        //                edges[i].ends[1].normal.x,
        //                edges[i].ends[1].normal.y,
        //                edges[i].ends[1].normal.z);
        //            f.Write("]\n\t\t \"P\" [ ");
        //            f.Write("{0} {1} {2} ", v1.x, v1.y, v1.z);
        //            f.Write("{0} {1} {2} ", v2.x, v2.y, v2.z);
        //            f.Write("{0} {1} {2} ", v3.x, v3.y, v3.z);
        //            f.Write("{0} {1} {2} ", v4.x, v4.y, v4.z);
        //            f.Write("]\n");
        //        }
        //    }
        //    f.Write("WorldEnd\n");
        //    f.Close();

        //    for (i = 0; i < vertices.getCount(); i++)
        //    {
        //        vertices[i].tex.x += 0.5;
        //        vertices[i].tex.y += 0.5;
        //    }
        //}

        #endregion

        //m is the vertex for which ring is the 1-ring
        //start is the first vertex in the ordered ring
        public UnicityList<Vertex> OrderRing(Vertex m, Vertex start, UnicityList<Vertex> ring)
        {
            int i, j;
            bool found;
            int miss = 0;
            Vertex v = start;
            UnicityList<Vertex> r = new UnicityList<Vertex>();
            ring.Remove(start);
            r.Add(start);

            while (ring.Count > 0)
            {
                found = false;
                for (i = 0; i < m.Faces.Count; i++)
                {
                    for (j = 0; j < m.Faces[i].Edges.Count; j++)
                    {
                        Edge e = m.Faces[i].Edges[j];
                        if ((e.Vertices[0] == v) && (ring.Contains(e.Vertices[1])))
                        {
                            v = e.Vertices[1];
                            r.Add(v);
                            ring.Remove(v);
                            found = true;
                            break;
                        }
                        else if ((e.Vertices[1] == v) && (ring.Contains(e.Vertices[0])))
                        {
                            v = e.Vertices[0];
                            r.Add(v);
                            ring.Remove(v);
                            found = true;
                            break;
                        }
                    }
                    if (found) break;
                }
                //no edge found, this could be a problem
                if (!found)
                    miss++;
                //the ring of vertices is probably broken
                // by a boundary, we should use reflection
                // to create virtual interpolation points.
                // instead we just look for the other edge of the boundary
                //it is not sufficient that this is a boundary point,
                // it must also only have one face in common with the
                // vertex for which the ring has been computed
                if (miss > 5)
                {
                    miss = 0;
                    for (i = 0; i < m.Edges.Count; i++)
                    {
                        if (m.Edges[i].IsBoundary)
                        {
                            if (m.Edges[i].Vertices[0] == m)
                                v = m.Edges[i].Vertices[1];
                            else
                                v = m.Edges[i].Vertices[0];
                            if (ring.Contains(v))
                            {
                                r.Add(v);
                                ring.Remove(v);
                                found = true;
                                break;
                            }
                        }
                    }
                }
            }
            return r;
        }

        //Compute the genus of the mesh (assuming all vertices and edges are used)
        //And assuming there are no boundaries
        //see Computer Aided Geometric Design (Farin) pg. 395
        public int SimpleGenus()
        {
            //Euler-Poincare formula
            return (vertexCount() - edgeCount() + faceCount() - 2) / 2;
        }

        public int faceCount()
        {
            return Faces.Count;
        }

        public int edgeCount()
        {
            return Edges.Count;
        }

        public int vertexCount()
        {
            return Vertices.Count;
        }

        public void RemoveFace(Face f)
        {
            int i;
            for (i = 0; i < 3; i++)
            {
                f.Vertices[i].Faces.Remove(f);
                f.Edges[i].Faces.Remove(f);
            }
            Faces.Remove(f);
        }

        public void RemoveFace(int j)
        {
            int i;
            Face f = Faces[j];
            for (i = 0; i < 3; i++)
            {
                f.Vertices[i].Faces.Remove(f);
                f.Edges[i].Faces.Remove(f);
            }
            Faces.Remove(j);
        }

        public int RemoveOrphanedEdges()
        {
            int i, n;
            Edge e;
            i = 0;
            n = 0;
            while (i < Edges.Count)
            {
                if (Edges[i].Faces.Count == 0)
                {
                    Edges[i].Vertices[0].Edges.Remove(Edges[i]);
                    Edges[i].Vertices[1].Edges.Remove(Edges[i]);
                    e = Edges[i];
                    Edges.Remove(i);
                    n++;
                }
                else
                {
                    i++;
                }
            }

            return n;
        }

        public int RemoveOrphanedVertices()
        {
            int i, n;
            Vertex v;
            i = 0;
            n = 0;
            while (i < Vertices.Count)
            {
                if ((Vertices[i].Edges.Count == 0) && (Vertices[i].Faces.Count == 0))
                {
                    v = Vertices[i];
                    Vertices.Remove(i);
                    n++;
                }
                else
                {
                    i++;
                }
            }

            return n;
        }

        private void CheckConnected(Edge e)
        {
            int i;
            Vertex v;

            //if (e==null) cerr << "Broken graph: edge" << endl;
            if (e.label == 1) return;
            e.label = 1;
            v = e.Vertices[0];
            v.label = 1;
            for (i = 0; i < v.Edges.Count; i++)
            {
                CheckConnected(v.Edges[i]);
            }
            v = e.Vertices[1];
            v.label = 1;
            for (i = 0; i < v.Edges.Count; i++)
            {
                CheckConnected(v.Edges[i]);
            }
        }

        private void CheckConnected(Face f)
        {
            int i;

            //if (f==null) cerr << "Broken graph: face" << endl;
            if (f == null) return;
            if (f.label == 1) return;
            f.label = 1;

            for (i = 0; i < 3; i++)
            {
                CheckConnected(f.Edges[i].Faces[0]);
                CheckConnected(f.Edges[i].Faces[1]);
            }
        }

    }

}

